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In  relerence  1 ol  the  report,  special  splittings  of  the  matrix  of  a sparse 
linear  system  were  proposed.  Ilvse  splittings  c.m  he  combined  with  the 
coniu-'ite  jtradient:'-  metltod,  which  results  in  von'  efficient  iterative 
so  In! ion  methods  when  the  matrix  is  a symmetric  M matrix.  The  research  report 
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(a)  a systematic  overview  of  possible  splittings  for  several  types  of 
matrices  and  eigenvalue  ini  omit  ion. 

( 

(bj  applications  to  nonsymnetric  linear  systems. 

(c)  solution  of  generalized  eigenvalue  problems  for  sparse  linear 
sys terns . 

\ 

The  research  on  (a)  lias  resulted  in  more  insight  on  how  to  choose  the 
appropriate  preconditioning  for  matrices  of  a special  structure.  Moreover, 
splittings  have  been  proposed  for  other  problems,  e.g.  30-problems  and 
periodic  boundary  condition  problems. 

The  research  on  (b)  has  not  yet  led  to  satisfactory  iterative  solution 
methods,  though  for  special  problems  an  always  converging  method  has  been 
deve loped . 

The  eigenvalue  problems  mentioned  under  (c)  have  been  solved  satisfactorily 
with  Lanczos-type  algorithms.  These  problems  were  met  in  the  investigation 
of  convergence  properties  of  the  methods  mentioned  under  (a)  and  (hi. 


. 
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Summary . In  ref.  I special  splittings  of  the  matrix  of  a sparse  linear 
system  were  proposed.  These  splittings  can  be  combined  with  the  conjugate 
gradients  method,  which  results  in  very  efficient  iterative  solution  methods 
when  the  matrix  is  a symmetric  M-matrix  l). 

The  research  reported  in  this  Final  Technical  Report  has  been  focused  on 
three  major  subjects: 

(a)  a systematic  overview  of  possible  splittings  for  several  types  of 
matrices  and  eigenvalue  information 

(b)  applications  to  nonsymmetric  linear  systems 

(c)  solution  of  generalized  eigenvalue  problems  for  sparse  linear 
systems . 

The  research  on  (a)  has  resulted  in  more  insight  on  how  to  choose  the 
appropriate  preconditioning  for  matrices  of  a special  structure.  Moreover, 
splittings  have  been  proposed  for  other  problems,  e.g.  3D-problems  and 
periodic  boundary  condition  problems. 

The  research  on  (b)  has  not  yet  led  to  satisfactory  iterative  solution 
methods,  though  for  special  problems  an  always  converging  method  has  been 
developed. 

The  eigenvalue  problems  mentioned  under  (c)  have  been  solved  satisfactorily 
with  l.anczos-type  algorithms.  These  problems  were  met  in  the  investigation 
of  convergence  properties  of  the  methods  mentioned  under  (a)  and  (b). 


')  A matrix  A “(a..)  is  an  M-matrix  if  a..v0  for  ifi,  A is  nonsingular  and 
-\  » .1  » .1 
A '0. 


Table  of  Contents 


Summary 

I.  Splitting  techniques  for  several  types  of 
symmetric  matrices 

II.  Iterative  methods  for  nonsymnetric  systems 

III.  Generalized  Eigenvalue  problems 

Literature  Cited 

Appendix  A 
Appendix  H 
Appendix  C 


Splitting  techniques  for  several  types  of  symmetric  matrices 


In  ref.  I the  idea  of  incomplete  decompositions  LU  for  arbitrary 
M-mat rices  A was  introduced.  For  symmetric  matrices  A those 
decompositions  were  used  as  preconditionings  for  the  conjugate  gradient 
algorithm  to  solve  the  linear  system  Ax»b.  In  the  examples  given  in 
ref.  I matrices  were  considered  arising  from  5-point  discretisation  of 
a self-adjoint  elliptic  partial  differential  equation  on  a rectangular 
region. 

During  the  period  of  the  Grant,  attention  has  been  given  to  a more 
systematic  investigation  of  possible  preconditionings.  From  the  information 
that  came  out  this  investigation  also  precodit ionings  could  be  proposed  for 
other  types  of  matrices.  We  mention  here  problems  that  come  from  p.d.e.'s 
with  periodic  boundary  conditions,  3D-problems  and  problems  where  the 
resulting  matrix  has  an  arbitrary  non-zero  structure.  Also  research 
has  heen  done  for  problems  where  the  matrix  is  not  an  M-matrix. 

This  research  has  been  reported  briefly  in  ref.  2 and  extensively  in 
Appendix  A.  As  a conclusion  we  mention  that,  with  respect  to  either 
economy  of  space  or  the  amounts  of  computational  work,  optimal  choices  for 
a preconditioning  can  be  made  for  problems  that  come  from  discretisations 
of  elliptic  self-adjoint  p.d.e.'s. 

The  eigenvalue  information,  reported  in  Appendix  A,  demonstrates  the 
effects  of  constructing  a more  or  less  incomplete  decomposition  for 
use  as  a preconditioning. 


II.  Iterative  methods  for  nonsymmetric  systems 


The  algorithms  mentioned  in  section  I work  all  fine  for  symmetric 
positive  definite  linear  systems  but  they  fail  to  work  for  nonsymmetric 
systems  since  they  are  based  on  the  conjugate  gradients  method. 

Variants  of  the  conjugate  gradients  algorithm,  that  converge  also  for 
certain  classes  of  nonsyimnetric  systems,  are  proposed  by  Concus  and 
Colub  [3],  Widlund  [4]  and  others.  They  share  the  disadvantage  that  a 
symmetric  system  has  to  be  solved  in  each  iteration  step.  Thus  in 
general  these  methods  can  not  compete  with  the  original  cg-method  for 
symmetric  matrices. 

It  has  been  investigated  during  the  Grant-period  whether  incomplete 
LU-decompositions  could  be  used  to  speed  up  the  convergence  of  some 
well-chosen  method.  Two  different  ways  have  been  followed,  each  of  which 
seems  to  have  promises  but  research  is  only  in  its  very  first  stage. 

The  first  way  was  to  base  iterative  methods  on  the  following  splitting 
of  a given  matrix  A: 

A = a(A  + AT)  + ((1  - a)A  - aAT) 

= aM  + N 

T T 

instead  of  A=j(A+A  )+j(A-A  ) which  is  more  usual. 

. T 

In  Appendix  B formulas  are  given  that  yield  a converging  method  if  A+A 

is  symmetric  and  positive  definite.  It  is  shown  in  Appendix  B that  the 

T 

replacement  of  M by  its  incomplete  LL  -decomposition  in  general  results 
in  faster  convergence. 

The  other  way  was  to  follow  ideas  expressed  by  Manteuffel  [5]  based  on 
Chebychev  acceleration  of  iterative  methods  for  non-linear  systems.  It 
was  recognized  that  the  use  of  incomplete  LU-decompositions  of  A could  be 
used  as  a preconditioning  and  they  gave  more  efficient  solution  methods 
than  for  instance  the  use  of  a Fast-Poisson-Solver  as  a preconditioning. 


III.  Generalized  Eigenvalue  Problems 


1 


In  order  to  study  the  convergence  properties  of  the  iterative  methods 
mentioned  in  I and  II  it  was  necessary  to  inspect  the  eigenvalue 
distribution  of  the  respective  iteration  matrices. 

The  Lanczos-algorithm  as  proposed  by  Paige  [6]  has  been  chosen  to  compute 
a large  number  of  eigenvalues.  In  ref.  7 the  numerical  experiments  are 
described  for  calculating  eigenvalues  of  BA  where  B and  A are  both 
symmetric  and  positive  definite. 

Since  also  eigenvalues  were  required  of  matrices  BC,  where  B is 

T 

symmetric  positive  definite  and  O-C  (antisymmetric),  a variant  of  the 
Lanczos-algorithm  is  proposed  in  Appendix  C. 

The  often  very  tedious  examination  of  the  numerical  results,  as  usually 
required  when  using  Lanczos-methods , has  been  largely  overcome  by  an 
automatic  selection  procedure.  This  research  is  also  described  in 
Appendix  C. 


Literature  Cited 


1.  Meijerink,  J.A.  and  van  der  Vorst,  H.A.,  An  iterative  solution 

method  for  linear  systems  of  which  the  coefficient  matrix  is 
a symmetric  M-matrix, 

Math,  of  Comp.,  Vol.  31  (1977),  pp.  148-162. 

2.  van  der  Vorst,  H.A.,  Approximate  decompositions  of  sparse  matrices. 

Proceedings  of  the  1977  Army  Numerical  Analysis  and  Computers 
Conference,  ARO  Report  77-3  (1977). 

3.  Concus,  P.  and  Golub,  G.H.,  A generalized  Conjugate  Gradient  Method 

for  Nonsymmetric  Systems  of  Linear  Equations, 

Proc.  2nd  Int.  Symp.  on  Computing  Methods  in  Applied  Sciences 
and  Engineering,  IRIA,  Paris  (1975). 

4.  Widlund,  0.,  A Lanczos  Method  for  a Class  of  Non-Hermit ian  Systems 

of  Linear  Equations, 

Courant  Institute,  New  York  University  (1976). 

5.  Manteuffel,  T.A.,  Adaptive  Procedure  for  Estimating  Parameters  for 

the  Nonsymmetric  Tchebychev  Iteration, 

Report  SAND  77-8239,  Sandia  Laboratories,  Livermore  (1977). 

6.  Paige,  C.C. , Computational  Variants  of  the  Lanczos  Method  for  the 

Eigenproblem, 

J.  Inst.  Math.  Applies.  10  (1972),  pp.  373-381. 

7.  van  Kats,  J.M.  and  van  der  Vorst,  H.A.,  Numerical  results  of  the 

Paige-style  Lanczos  method  for  the  computation  of  extreme 
eigenvalues  of  large  sparse  matrices, 

TR-3,  Academic  Computer  Center  Utrecht  (1977). 


Guide  lines  for  the  usage  of  incomplete  decompositions 
in  solving  sets  of  linear  equations  as  occur  in 
practical  problems. 

by 

J . A.  Meijerink  + and  H.A.  van  der  Vorst. 


+■ 


Koninkl i jke  Shell  Exploratie  en  Produktie  Laboratorium 
Rijswijk,  The  Netherlands. 


IT 


CONTENTS 

Page 

III 

Summary 

1.  Introduction  1 

2.  Incomplete  decompositions  for  symmetric  M-matrices  2 

. .1  rive-point  discretisation  of  elliptic  p.l.e.'s  in  2-D  2 

2.1.1  Diagorwil  scaling  3 

2.1.2  ICCG(1,1)  and  SSOR  OJ=l)  3 

2.1.3  ICCG(1,2)  4 

2.1.4  ICCG(1,3)  4 

2.1. 5 ICCG( 2,4)  5 

2.1.6  Some  other  decompositions  6 

2.2  Five-point  discretisation  for  problems  with  periodic 

boundary  conditions  6 

2.2.1  ICCG( 1,1)  6 

2.2.2  ICCG( ] ,2)  7 

2.2.3  ICCG(l,l)~per iodic  8 

2.3  Seven-point  discretisation  of  elliptic  p.d.e.'s  in  3-D  9 

2.3.1  ICCG(1,1  jD  9 

2.3.2  Other  decompositions  for  3-D  9 

2.4  K-matrices  aidsing  from  five-point  discretisations  on 

irregular  regions.  10 

2.5  M-ma  trices  with  an  irregular  non-zero  structure  11 

3.  Algorithms  for  symmetric  positive  definite  matrices  11 

4.  Algorithms  for  non-symmetric  positive  definite  matrices  12 

5.  Numerical  experiments  12 

References  14 


Tables  1 and  2 
Figures  1-21 


1 


III 


SUMMARY 


This  report  presents  incomplete  decompositions  for  various  types  of 
matrices  as  occur  in  the  implicit  discretisation  of  practical  problems. 

A review  is  given  of  methods  for  the  usual  five-point  discretisation  of 
a self-adjoint  elliptic  second-order  partial  differential  equation  in  two 
dimensions  on  a square.  The  matrices  which  occur  in  this  type  of  problem 
are  symmetric  M-matrices  of  very  regular  structure.  The  convergence  behaviour 
of  the  different  decompositions  for  this  case  is  demonstrated  by 
numerical  experiments.  The  report  also  gives  decompositions  for  the 
following  type  of  matrices: 

(i)  Symmetric  M-matrices  of  a different  structure 

(ii)  Symmetric  positive  definite  matrices 

(iii)  Non-symmetric  positive  definite  matrices 
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GUIPK  LINKS  l "OR  'me  USAGE  OK  INCOMTUT'K  nKOOMl'OSIT IONS  IN  SOLVING  SETS 
OK  LINEAR  EQUATIONS  AS  OC'CUK  IN  PRACTICAL  PROBLEMS 


l.  INTROPUCTION 

In  Rot.  1 the  idea  of  incomplete  decomposition  U.J  for  arbitrary  M-matrices 
A war.  introduced . Kor  symmetric  nutrices  A tltose  dtvompositions  were  used  as 
precondit  ioniiy.s  toi’  the  conjugate  gradient  algorithm  to  solve  the  linear 
equation  Ax-b.  In  the  examples  given  in  Ref.  1 nutrices  were  considered  arising 
trom  S-point  d i sore t iso  l ion  oi  a self-adjoint  elJipiic  j>irt  ial  differential 
equation  on  a rectangular  region.  Only  ti\V>  different  incomplete  decompositions 
were  demons t rn  t <\1 . 

In  this  report  we  present  a more  systemit ic  review  of  the  possible  incomplete 
do.-omu':  i t ions  tor  that  problem  (section  1 . 1 1 . In  addition,  we  shill  discuss 
incomplete  decompositions  tor  other  types  of  matrices,  e.g.  M-matrices 
arising  from  problems  with  periodic  boundary  conditions  (section  3.?) 
and  M-iiut rices  with  a nor'c  arbitrary  structure  (sections  3.4  and  3.S).  for 
this  purpose  we  need  an  extension  of  the  definit ion  of  an  incomplete 
decomposition  g.iven  in  the  prwt  ot  theorem  2.3  in  Ret.  1.  In  the  process 
defined  there  some  off-diagonal  elements  were  omitted  after  each  elimination 
step.  It,  instead  ot  omitting  some  off-diagonal  elements,  we  replace  them 
by  neg.it  ive  elements  which  are  smaller  in  absolute  value,  or  if  diagonal 
elements  are  replaced  by  larger  ones,  the  construction  process  does  not 
break  down  and  the  resulting  incomplete  decomposition  'UK  defines  a regular 
splitting  of  A.  This  new  process  describes  the  extended  concept  oi  incomplete 
decompositions.  A specific  example  of  this  process  is  the  following:  In  the 
k^1  step  of  the  Gaussian  elimination  pvvess  elements  are  eliminated  with 
the  k**’  row.  This  nuy  cause  three  effects: 

i)  :’,ero  off-d iagotvil  elements  turn  to  negative  non-zero  values 

ii)  non- zero  off-diogoml  elements  hvomr  snviller  (oltliough  larger  in 
absolute  value) 

i i i )d  iagou.il  elements,  Kwnne  smaller  (but  amiin  iv>sitivo). 

Thus  omitting  to  carry  out  the  elimiiwtion  correct  ions  tor  some  of  the 

elements  ot  the  nutrix  results  in  an  incomplete  deoompos i t ion.  Examples  of 

this  type  of  deoompos. i t ion  are  given  in  soot  ions  2.1.2,  2.1.3,  2.4  and  2.S. 

2 3 

Other  opproximite  too  tori  sot  ions  arc  discussed  by  Stone  , Itupont  ot  all  and 


l( 

dust  aft  son  . They  all  mttwluce  one  in-  more  jvirameters  into  the  decompor.  1 1 ion 
1'iu'i‘s-;  to  accelerate  convergence.  The  pivperty  ot  regular  splitting  is  lost 
in  uv'st  ot  their  examples. 

hershiw''  and  Mant  out  t el1'  [u\'v ide  extensions  tor  positive  delinite  na trices . 

We  skill  describe  these  hrictly  ami  present  another  one  in  sect  ion  3. 
Incomplete  deoompos i t ions  tor  p.d.e.'s  in  three-dimensions  are  treated  in 
sivt  ion  '3.3.  In  sect  ion  li  algorithms  tor  non-symmetric  nvitr ices  ore  described. 
In  section  S convergence  results  as  well  as  'eigenvalue'  Lnfomat  ion  on  the 
deeompos i t ions  described  in  sect  ion  .’.1  are  given  for  some  speeijic  examples. 


;.  iNtoMi'i.nr  PixvMiw.moNS  tor  synnitru'  m-matruts 


.’.1  Tiyej|vdjU  discret  isat  ion  ot  elliptic-  pirlial  d i 1 1 erent  ia  1 ogiut  ions  in 
two-d  iniens  ions 

The  I inear  equal  ions  in  this  sect  ion  arise  1 1 v 'in  t ive-polnt  discrete  apptwxi- 
nvit  ion  to  the  soi-ond-omler  selt  adjoint  elliptic  jvirt  iul  ditterential  txjn.it  ion: 

£ At x,y)  u(x,y)  - IKx,y)  ^7  u(x,v)  ► dlx,y)  ulx,y)  = lKx,y) 

(2.1.1) 

with  A(x,y),  lKx,y)-»0,  d(x,y)>l)  ami  (x,y)  r K,  where  K is  a rectangular 
region.  A l oi ig  the  loumlury  iSR  ot  K the  h'undary  condition 

alx.y)  ulx,y)  * U(x,y)  y^u(x,y)  = >(x,y) 

tk'lds,  with  n(\,y),  tUx,y)'0  amt  a(x,y)  ♦ t^lx,y)'h  ami  wliere  is  the 
outw.uvl  derivat  ive  perpend icular  to  tSK . Tin-  s.tructure  ot  the  resulting 
synnietrie  M-matrix  A ot  oivler  N is  shown  in  Tig.  1.  The  elements  ot 
tin’  diagouil  ot  A are  denoted  by  a . , t lose  ot  the  t irst  upper  diugoivil  by  b- 
and  those  ot  the  m-tli  upper  diagouil  by  c . , where  i is  the  index  ot  the  iv>w 
ot  A in  which  the  respective  elements  cxvur,  and  m is  the  hilt  lundwidth 


ot  t lie  matrix.  Tor  the  di'rivat  ion  ol  such  l inear  systems  see  Ret  . 7. 


r w 
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2 . 1 . 1 _U  ing 

1'he  simplest  allowed  ctioiee  for  K is  the  diagoivil  of  A.  The  resulting 
conjugate  gradient  method  is  the  same  as  the  c.g.  method  applied  on  the 
matrix  scaled  by  its  diagonal.  This  seal ing  is  in  some  respect  optimal,  since 

it  minimises  approximately  the  condition  number  of  h among  all  diagonal 

. . 3 

scalings  . 

If  the  equation  is  scaled  in  advance,  the  number  ot  multiplications  is 
10  N per  iteration.  It  not,  this  number  will  be  11  N. 


2.1..’  ICCG_ ( 1 ^ n _a nd_SoOR_ (w = 1 ) 


T 


Here  the  matrix  K is  chosen  so  that  its  decomposi lion  factor  1.  has  the 
same  sparsity  pattern  as  the  upper  triangular  port  of  A.  This  decomposition 


lus  been  considered  by  many  authors 


12  t 4 0 


= U 


D. 


It  is  convenient  to  write  this  decomposition  as  k.  , - l . . , o . i..  . , 

i|>  l » 1 1 * l 1 > L i i 1 

where  h.  . is  an  upper  triangular  matrix  and  D,  a di.agoiul  matrix  equal 

-L  > L ,p  L y L 

to  the  inverse  of  the  main  diagonal  of  L.  . . In  common  with  the  elements  of 


T 


M' 


A,  those  of  L . are  denoted  by  a,  b.  and  e.  and  tlx'ise  of  D.  by  d. 

i y l 1 -L  1.  1 ’ 1 

(see  fig.  .').  The  following  relations  axv  easily  verified: 


a.  = d.  ' = a.  -b-‘,  d.  . -c.‘  d- 
l l a l-l  t- l i~m  l-m 


b.  = b.  and  c-  = c. 
it  li 


for  i = l ,2, . . . ,N 


In  these  relations  the  non-dot  Lnod  elements  .ire  zero.  Only  extra  s toixage  for 
the  elements  d-  is  required.  The  resulting  hybrid  conjugate  gradient  method 
is  called  I0CG(1,2).  The  indices  are  used  to  indicate  that  there  is  one 
non- zero  diagonal  next  to  the  main  diagonal  and  one  non-zero  diagonal  at  the 
outward  side  of  the  hand.  This  is  1CCG(0)  of  Kef.  1.  The  number  of 
multiplications  for  this  me that  is  16  N multiplications  per  iteration. 

The  SSOR(ui-l)  decomposition  arises  if  all  Gaussian  elimination  corrections 
are  neglected.  Thus  SS0R(w=l ) is  an  example  ot  the  extended  class  of 
incomplete  deconipos i t ions  mentioned  in  section  1.  The  nimiber  of  multi- 
plications remains  the  same  as  for  IC0G( 1,1),  but  one  array  of  storage 
his  been  saved.  For  the  use  of  SSOR  as  a precoixlit ioniug  technique  see 

run. 


1 


]vAAj(l12) 

lh'“  llv*lrAX  4,1  4,1  4,1  4,1  ot  tf,e  I'i'evioua  section  is  .4  matrix  equal 

li’  A,  except  tor  two  diagonals  adjacent  to  the  outermost  two  diagonals, 

• is  indicated  by  the  dotted  lines  in  Fig.  .1.  Hy  including  non-zero  entries 
on  ttK'se  lines  in  L and  L , we  e\|>eot  to  Lmprvn'e  the  approximate  decompos i t ion . 
Hus  approximation  will  he  written  as 

h.  h n 1 * 

l-’  1,2  1,2  4,.> 


where  l'^.,  is  the  d iagoiwl  matrix  equal  to  the  inverse  of  the  mtin 


diagonal  of  the  elements  ot  l.J  , are  denoted  by  a.,  b.  , . and  r. 

and  those  of  D1j2  by  ch  (see  Fig.  U).*  lire  elements  a.,  fV,  6.,  ,V  and  e. 
vMu  K‘  computed  ixvursively  tiv'in 


i ‘ i ‘l  i-l  ' i-tn+1  4i-m+l  ' i-m.  4-m 


h.  b.  - o.  d.  e 

1 1 i-rtt+1  i-m+ I i-m+1 


for  i l ,?,... ,N 


'i-l  4-1  4-1 


Ihe  ix'n-det  ined  elements  ai'e  all  .-era. 

Storage  is  required  tor  three  arrays  ot  length  N.  lire  number  of  multipli- 
cations necessary  for  each  iteration  step  ot  the  resulting  hybrid  conjugate 
gradient  met  tied  ICCG(1,2)  is  equal  to  IS  N. 

In  otx|er  to  save  computer  storage  (one  array)  the  relatively  snail 
wnnvction  on  b.  can  be  omitted.  'thus  b.  = b..  This  is  another  example  ot  tire 
extended  class  ot  incomplete  d.vxrmpos it  ions  and  will  be  denote!  by  li\V,  vl*."' 

Ihe  matrix  , it:  equal  to  A,  except  tor  the  two  dottel  diagonals  as 
indicated  m Fig.  b.  'these  non-zero  elements  can  lv  eliminated  by  introducing 
an  extra  diagorvtl  in  1.  (see  t ig.  ft).  lh«*  elements  on  these  d iagouils  will 
be  denotixl  by  t..  Ill  is  incomplete  dexompos  i t ion  will  U>  denote!  by 

K,  , s l,  , p 1 1 

M 1,3* ’.,3  4,3 


a, 


I llu*  eleircentu  of  l'.  and  L.  ' . can  l •« • computed  livm: 

1 l , -S 


I 


f> 


2.1.  t > Uome_other  deeompos  it  ions 

Proceeding  in  this  tanner,  we  obtain  a sequence  of  incomplete  decompositions 
Kj  0,  Kg  ^l(,  , ,,  etc.  From  this  we  see  that  the  number  of 

non-zero  diagonals  grows  rapidly  and  thus  the  amount  of  work.  However’, 
for  most  practical  problems  the  two  diagonals  next  to  those  of  the 
previous  decomposition  cause  the  major  improvement.  In  this  way  K(3,5) 
h(4,b)  etc.,  together  with  the  corresponding  ICCG  methods  ICCG(3,5), 

ICCG(4,6)  etc.,  are  developed. 

Up  to  row  we  have  always  added  complete  diagoruls  in  the  decomposition 
process.  The  diagonals,  however,  contain  their  non- zero  entries  only  in  a 
blockwise  structure  (see  Fig.  9).  In  part icular,  this  implies  that  in  our 
terminology  a complete  choleski  factorisation  is  equivalent  to  "incomplete 
decomposition"  with  2m- 2 extra  "diagonals" . 

2.2  Five-point  discretisation  for  problems  with  periodic  boundary  conditions 

The  linear  equations  in  this  section  arise  in  the  same  way  as  the  linear 
equations  in  section  2.1,  except  tint  a periodic  boundary  condition  holds  in 
at  least  one  direction.  In  the  examples  we  have  restricted  ourselves  to  a 
periodic  boundary'  condition  in  the  x-direction.  Tins  boundary  condition  gives 
rise  to  additional  elements  in  the  matrix  A,  as  indicated  in  Fig.  10.  These 
£xtra  elements  are  denoted  by  6^.  Since  i is  the  row  index,  only  3^ , 

^2m+l’’*'  are  Again  A is  an  M-mvtrix. 

2.2.1  ICCGQ^l) 

In  common  with  the  non-periodic  case  in  2.1.1  air  incomplete  decomposition 
can  be  constructed  in  cases  where  the  upper  triangular  factor  lus  the 
same  sparsity  structure  as  the  upper  triangular  part  of  A.  This  decompo- 
sition is  written  as 

Kl,l  = 4,1  Dl,l  4,1 


The  elements  of  L,  , arid  D . can  be  calculated  from 
L ? L 1 » I 

- - -1  ~ 9 ~ 9 - -9 

a.  = d.  = a.  -b.  . d.  . -8-  .-id-  -c  • d. 

ii  l i-I  l-l  i-m+1  i-m+1  l-m  l-m 


b.  = b.  c.  = c.  8-  = 8 • 

ii  ii  ii 


i=1.2,...N 


Non-defined  elements  are  zero.  Note  that  the  term  i 0 only 

if  i=m,  2m,  3m,....  The  resulting  ICCGC1,1)  process  again  takes 
approx imately  16N  multiplications  per  iteration. 

2.2.2  IC0G(lj2) 

The  matrix  K.  . of  2.2.1  has  elements  as  indicated  in  Fig.  11.  To  annihilate 

X J 1 rp 

these  elements  in  L , non- zero  elements  are  required  in  these  places.  These 

elements  are  denoted  as  shown  in  Fig.  12  by  a^,  to,  c^,  £h,  e^,  y^, 

and  the  elements  of  D „ by  d . . They'  can  be  calculated  from 

i , z 1 

~ ~ -1  '?  ~ 2 ~ 2 
a.  = d-  = a.  -b.  . d.  . -e.  d.  ..  -c.  d. 
ii  l i-I  i-l  i-m+1  i-m+1  l-m  l-m 


- 6.  , d.  . 
i-l  i-l 


(only  for  i=m+l,  2m+l,  3m+l, — ) 


~ 2 ~ ~ 2 

- y.  d.  -8-  ...  d.  , (only  for  i=m,  2m ) 

i-m+2  i-m+2  i-m+1  i-m+1  J ’ 


b.  = b.  -c.  . d.  ..e.  . 

l l i-m+1  i-m+1  i-m+1 


ci  = ci 


e-  = -c.  . d.  , b.  . 
l i-l  i-l  i-l 


6.  = -c.  „ d.  y.  -c.  .nd.  ,i8-  i=m,  2m — 

l i-m+2  i-m+2  i-m+2  i-m+1  i-m+1  i-m+1 


i=l,2,3...N 


y.  = -8-  , d.  . b.  i 
l i-l  i-l  i-l 


i=2,  m+2,  2m+2, 


8 • = 8 • -c . . d . . 6 . . 
l l i-l  i-l  i-l 


i=l,  m+1,  2m+l. 


Note  tlvit  b , b,,  , — and  e.,  e are  non-zero. 

m’  2m’  1’  m+1’ 

The  resulting  ICCG  process  takes  1 BN  multiplications  per  iteration  and 
requires  roughly  three  arrays  of  length  N for  the  incomplete  decomposition. 


8 


L 


2.2.3 


The  incomplete  decomposition  of  2.2.1  does  not  hive  a periodic  structure. 
To  obtain  a K with  a periodic  structure  we  write 


K -CL  + D _1 ) D (l7  + D"1) 

p - - - 


P P P P P 


The  periodic  structure  of  the  matrix  L is  given  in  Fig.  13.  D is  a 

p B & p 

diagonal  matrix.  The  elements  of  L and  D tave  to  satisfy: 

P P 


b.  = b. 
l i 

c.  = c. 

11  11 

for 

i=l. 

2,-N 

(2. 2. 3.1) 

= a:1 
1 1 

' 2 ' ~ ' 

= a . -b . . d . , -c . d • 
i l-l  l-l  l-m  l-m 

for 

i=2. 

3, — , m,  m+2,  m+3, — 

,2m,2m+2. 

. . . ,N 

(2.2.3 .2) 

= 5:1 

1 1 

= a.  -0?  d.  . -c.‘  d. 
i i i+m-1  l-m  l-m 

for 

i=l. 

m+1,  2m+l, — , N-m+1 

(2. 2. 3. 3) 

The  d.  cannot  be  calculated  straightforwardly,  since  in  the  second  formula 
di+.,_l  is  present.  We  can  calculate  them  by  substituting  (2.2. 3.3)  into  (2. 2.3.3) 
for  the  next  i,  and  continuing  in  this  way,  we  find  quadratic  equations  for 


the  d^.  For  we  choose  the  largest  root,  since  this  choice  results  in 

smaller  elements  b . d . ' c . , in  the  error  matrix  K -A.  We  now  give  the 

l l-l  r-m-1  p c’ 


derivation  for  the  formula  for  d . The  d.  can  be  computed  in  a similar'' 

m km  1 


way . We  rewrite  ( 2 . 2 . 3 . 3 ) as 


d. 1 = w.  - v . d 
1 1 l m 


(2. 2. 3. 4) 


and  2. 2. 3. 2 as 


d7  = w.  - v.  d.  , 

l li  l-l 


i=2, m 


(2.2.3.b) 


From  induction  it  follows  that: 

p.  + q.  d 
a m 


d.  = 

v . + s • 3 


for  i=l,  2,  — ,m 


i m 


Hie  coefficients  p^,  qj  and  s^  satisfy  p^  = 1 = 0 r^  = w^  s^  = -v^ 


pi+l  = ri  qi+l  = Gi  = wi*i  ri  P,-  s,ai  = s,--~  vmq 


■i+1  "i+1  *i  ’i+1  Fi  i + I ' wi+l  i i+lsi 

IT. is  leads  to  the  quadratic  equation  in  with  known  coeeficients  p^,  q^. 


r and  s : 
rr  m 


2. 3 Seven-point  discretisations  of  elliptic  p.d.e.'s  in  three  dimensions 

The  seven-point  discretisation  for  equation  (2.1.1)  in  three  dimensions 
leads  in  a similar  way  to  a matrix  with  seven  non-zero  diagonals.  The 
structure  of  this  matrix  A is  shown  in  Fig.  14.  The  elements  of  the  upper 
triangular'  part  of  A are  denoted  by  a^,  b^,  and  e^,  where  i is  counted 
rowwise.  If  n,  m,  k are  the  number  of  gridpoints  in  the  x,y,z  directions, 
respectively,  the  order  of  the  matrix  and  the  sizes  of  the  blocks  are 
nxmxk,  nxm  and  n. 


2.3.1  ICCGC 1,1,1) 

In  common  with  the  2-D  case  the  ICCG  (1,1,1)  factorisation  is  the  one 

where  the  upper  triangular  factor  has  the  same  non-zero  structure  as  the 

upper  triangular  part  of  A.  Again  this  decomposition  is  written  as 

T T 

Ki  i i = Lf  i ^ ^ ^ , where  is  an  upper  triangular  matrix 

and  D1  . . a diagonal  matrix  equal  to  the  inverse  of  the  main  diagonal  of 
L . . The  elements  of  L , , are  denoted  by  a.,  To • , c . and  e . and  the 

1,1,1  _ T.,1,1  J i’  l*  i l 

elements  of  ^ by  cL . 'Ihese  elements  are  given  by  the  recurrency 

relations: 


for  i=l,2, ,n  x m x k 


Non-defined  elements  should  be  replaced  by  zeros.  It  can  easily  be  seen 
that  for  major  problems  where  the  diagpnals  cannot  be  stored  all  together 
in  core,  the  d^  can  be  calculated  by  taking  successively  only  parts  of  the 
order  of  n x m in  core. 

The  resulting  hybrid  conjugate  gradient  method  requires  20  N 
multiplications  per  iteration. 
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2.3.2  Other  deeomro s i ^ io  ns  _ t or  _ J -D 


Hie  matrix  K.  . , = L.  . .P.  . . L.  . . , of  the  pivvious  section  is  a 
natrix  equal  to  A,  except  lot'  six  diagonals,  as  shown  in  Fig.  15)  as  dotted 
lines.  We  obtain  the  next  decomposition  by  including  non- zero  entries  on 
these  lines  in  L and  L1 . The  elements  of  Ll  are  denoted  by  a^,  In  , , e^ , 

f.,  g.  and  h.,  as  shewn  in  Fig.  16,  and  can  be  calculated  from: 

i 6i  l 


a . * U . * = a • - bt  . 1 . . 
i i i i-l  i-l 


N-n*l  l'i-n*l 


J. 

l-n  l-n 


i-nm*n  i-rw,*n 


-t  . 


i-im*l 


i-ra*l 


- 2 
1 i-nm 


a. 

i-nm 


b. 

i 


b.  -c.  , d.  , h. 

1 1-11*1  1-11*1  1-11*1. 


-(•  (J . f 

i-mn*l  i-mn*l  l-tmi*l 


1). 

— 0 * . d • % 
l-l  1-1 

1 

. 

_ 

c . 

= C . -e . 

1 

l i-mi 

= -f.  , d 

i-n*l 

f. 

= -e.  , d. 

i 

i-l  l- 

e. 

- tf  • 

i=l, . . . .ink. 


, v* . , h ■ , -«*  ■ <i  - c - 

i-ii*l  i-ii*l  i-n*l  i-n  i-n  l-n 


Six  arrays  of  length  N are  inquired  to  store  the  non- zero  diagonals  of  L . 
Hie  resulting  ICCG  metliod  ivquires  26N  multiplications  per  iteration. 
Unfortunately,  it  we  pi  weed  in  this  manner  the  number  of  diagpruls 
in  the  subsequent  decompositions  will  increase  rapidly.  For  instance,  the 
next  decomposition  has  12  non-zero  diagonals  in  its  upper  triangular  part, 
lire  resulting  ICCG  method  takes  36  N multipl ications  per  iteration. 


3.U  M-nvatrices  arising  from  five-point  discivt  i sat  ions  on  irregular  regions 

So  far  we  tvive  only  considered  discretisations  on  square  regions.  We  are 
now  going  to  conrnont  on  regions  with  internal  boundaries  (no- flow 
hourkLiri.es>  or  differently  shaped  regions.  For  convenience  it  will  he 
assumed  th.it  the  region  consists  of  snvil  l squires. 

An  .internal  boundary  is  reflected  by  some  extaw  zeros  in  the  matrix,  but 
the  matrix  remains  a symmetric  M-nutrix,  thus  incomplete  decompositions 
can  be  constructed  as  before.  An  internal  Kx lidary  implies  that  there  is  no 
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direct  connection  (no  flow)  between  points  on  different  sides  of  the 
haundary.  This  property  is  preserved  in  each  of  the  above-mentioned  decompo- 
sitions.  This  is  in  contrast  with  Stone's  SIP  method  , where  the  use  of  the 
iteration  parameter  may  cause  a connection  through  a no-flow  boundary. 
Irregularly  shaped  regions  can  be  extended  in  an  obvious  way  to  square 
'regions  with  an  internal  boundary  at  the  point  of  the  original  real  boundary. 
The  linear  system  arising  from  this  extended  region  can  be  treated  as 
before,  bearing  in  mind  that  the  extended  parts  do  not  require 
computational  work. 


2. 5 M. matrices  with  an  irregular  non-zero  structure 

M-ma trices  with  an  irregular  non-zero  structure  arise,  for  instance, 

from  some  finite-element  methods  on  irregular  meshes^  and  pipeline  net- 

, 12 
rocks 

We  write  the  matrix  AasA  = L + D + U where  L,  U are  strictly  lower 
and  upper  triangular,  respectively,  and  D the  diagonal  of  A.  If  we  omit 
all  Gaussian  elimination  corrections  on  off-diagonal  elements  (see  section  1), 
then  the  incomplete  decomposition  is  given  by  = (L  + D^)  (D  + U). 

is  determined  by  the  relation  tliat  the  diagonal  of  Kq-A,  which  is 
equal  to  the  diagonal  of  Do  + diag(LDQ  \l)  -D,  is  zero. 

If  the  nutrix  is  symmetric,  this  decomposition  can  be  combined  with 
the  conjugate  gradient  method.  For  non-symmetric  matrices  see 
section  4. 


3.  ALGOR  TIT  iMS  FDR  SYMETRIC  POSITIVE  DEFINITE  MATRICES 


If  the  natrix  is  not  an  M-matrix,  the  construction  of  an  incomplete 
decomposition  nay  fail,  because  of  the  occurrence  of  non-positive  diagonal 
elements'  . Snvill  positive  diagonal  elements  are  also  undesirable,  because  of 
stability  problems . 

Three  different  strategies  which  seem  to  overcome  this  problem  have 
currently  been  proposed: 
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(.i)  If  a diagonal  element  of  less  than  a prescribed  positive  value 

T 

is  encountered  during  the  construction  of  the  incomplete  LL 
decomposition  then  some  already  computed  off-diagonal  elements 
in  the  corresponding  column  of  L are  set  to  zero. 

t i i ) The  diagonal  element  can  be  enlarged  if  necessary  by  adding  a 
sufficient  amount  to  the  original  element^. 

tiii)  We  can  also  add  a I to  the  matrix  . If  a is  large  enough,  the  signalled 
problems  will  not  occur . 

4.  AIGPKLTHHS  FOR  NON  SYMMETRIC  POSITIVE  DEFINITE  MATRICES 

If  the  matrix  is  non- symmetric,  then  an  incomplete  LU  decomposition  K 
can  be  constructed  in  a similar  way  as  described  previously  for  the 
symmetric  matrices.  Since  symmetry  and  positive-definiteness  are  both 
required  for  the  conjugate  gradient  algorithm,  the  CG  algorithm  can  be 
applied: 

T T-^  -1  T T~^  -1 

A K K A x = A K K lb 

This  algorithm  requires  twice  as  much  work  per  iteration  as  the 
corresponding  symmetric  case  and  the  upper  bound  for  the  number  of 
iterations  increases  by  a factor  of  two. 

5.  NUMTRICAL  FXFLKINENTS 

To  obtain  an  impression  of  the  convergence  behaviour  of  different 
incomplete  decompos it  ions , we  have  for  the  ICCG-methods  introduced  in 
section  2.1: 

(i)  compared  the  convergence  results 

(ii)  calculated  the  eigenvalue  distribution  of  the  preconditioned 
mi  trices  K ^A. 
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The  two  test  problems  were: 

i)  problem_l_1  The  five-point  discretisation  of  the  Fbisson  equation  Au=0 

over  Osxsl,  0;>ysl  with  boundary  conditions  ^ = 0 for  x = 0 and  x = 1, 

3u  dx 

^77  - 0 for  y = 1 and  u = 1 for  y = 0.  A uniform  rectangular  mesh  was  chosen, 
with  Ax  = 1/31  and  Ay  = 1/31,  which  resulted  in  a linear  system  of  992 
equations.  The  solution  of  this  equation  is  known  to  be  u(x,y)  = 1 and 
as  initial  starting  vector  for  the  iterative  schemes  a vector  was  chosen 
with  all  entries  random  between  0 and  1.  This  was  done  to  prevent  co- 
incidental fast  convergence. 

...  7 

n)  problem-^  This  problem  has  been  taken  from  Varga  . Equation  (2.1.1) 

holds  for  R,  where  R is  the  square  region  CKx,  y<2.1  as  shown  below 


On  the  houndary  of  R,  the  boundary  conditions  are  ~~  - 0.  Further, 
D(x,y)  = 0 over  R and  the  functions  A,  B and  C are  given  by 


Region 

A(x,y) 

B(x,y) 

C(x,y) 

1 

1.0 

1.0 

0.02 

2 

2.0 

2.0 

0.03 

3 

3.0 

3.0 

0.05 

A uniform  rectangular  mesh  was  chosen  with  0.05  mesh  spacing,  so  that  a 
system  of  1849  linear  equations  resulted.  The  solution  of  the  system  is  known 
to  be  u = 0.  A vector  similar  to  the  one  in  problem  I was  chosen  as  a 
starting  vector. 

In  Tables  1 and  2 the  convergence  results  are  listed.  In  both  tables  we 
can  conclude  from  the  last  column  that  ICCG(1,3)  is  almost  optimal. 

Since  the  convergence  behaviour  depends  on  the  eigenvalue  distribution, 
where  the  coni  it  ion  number  uni  clustering  play  an  important  role, 
a number  of  the  lowest  and  highest  eigenvalues  have  been  calculated. 
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The  eigenvalues  are  all  divided  by  A^^  for  the  matrix  of  problem  1, 

preconditioned  with  several  incomplete  decompositions  because  an 

upper  bound  for  the  convergence  of  the  conjugate  gradient  method  is 

given  by  (/c^D/Cn^l) , where  c = A /A  . . The  distribution  of  these 

max  min 

scaled  eigenvalues  has  been  plotted  in  Figs.  17-21. 
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Iterative  methods  for  a class  of  nonsymmetric 
systems  of  linear  equations,  based  on  splitting-off 
a symmetric  part. 


by 


H.A.  van  der  Vorst. 
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Iterative  methods  for  a class  of  non-symraetrlc  systems 
of  linear  equations,  based  on  splitting-off  a symmetric 
part . 


Henk  A.  van  der  Vorst 


1 . Introduction. 


In  recent  papers  Concus  & Golub  /^\ / and  Widlund  / 2/ 
discuss  conjugate  gradient  like  iterative  methods  for 
the  iterative  solution  of  real  non— symmetric  algebraic 
equations . 

A x = b (l  .1  ) 

These  methods  are  based  on  the  splitting  of  the  matrix 
A in  a symmetric  and  a skew-symmetric  part 

A - £ (A  + AT)  + A (A  — AT)  (1.2) 

T 

and  the  requirements  are  that  ^ (A  + A ) is  positive 
definite.  In  this  paper  a class  of  splittings  of  the 
matrix  A is  considered  and  the  influence  of  the  special 
choice  of  these  splittings  on  the  convergence  of  a 
simple  related  iterative  method  is  discussed. 

A combination  of  these  splittings  and  the  incomplete 
choleski  factorization,  described  by  Meijerink  & van  der 
Vorst  /4/,  has  been  treated  in  more  detail.  Simple 
numerical  experiments,  showing  the  various  effects  are 
demonstrated.  Methods  for  acceleration,  such  as  those 
based  on  Manteuffel's  ideas  / 3/,  are  not  considered  here, 
although  they  might  be  very  effective,  since  the  eigen- 
values of  the  occurring  iteration-matrices  are  located 
in  the  complex  plane  on  straight  lines  parallel  to  the 
Y— axis . 


-2- 


2.  Splitting-off  a symmetric  part. 


'.ttf  C ' | i<ol(vj 

of  raatrlcnn  that 


Throughout  this  paper  we  think  of  matrices  that 

arise  in  the  five-point  discretisation  of  elliptic  partial 

differential  equations,  that  have  essential,  i.«.  not 
4 *■  l ■ t v \ 

OMT  removable,  first  order  differential  terms.  Let  us 

for  simplicity  think  of  the  equation 

A U ■»  = 0 (p.l) 

(in  this  case  however,  the  first-order  term  is  easily 

removed,  if  ^ is  some  continuous  function). 

i’or  the  matrices  A,  arising  in  thin  way,  it  follows  that 
T 

A + A is  a symmetric  M-matrix,  and  thus  positive  definite 
(see  Varga  /r)/). 

Very  simple  splittings  that  take  advantage  of  this  are 
splitting  I A - f(A  + A r)  + 1 ( A A1’)  (?,?.) 


splitting  11 


m rn 

A - (A  + A ) - A1 


(2.3) 


Let  ua  now  consider  a simple  basic  iterative  method,  that 
arises  from  the  splitting  A » M - N i 


M x.  , - b + N x. 
i-f  1 i 

This  leads  for  splitting  1 and  11  reap,  to 


(2.4) 


(2.3) 


^ \ (2.6) 
for  the  solution  of  (l.l). 

If  the  respective  errorvectors  and  U(  are  defined 

by  _ 

\ *i  * * ‘ and  C;  . \i-  * (2.7) 

where  x is  the  solution  of  Ax  b,  then  we  have  the 


rein t ions 


V.  - (O.rt- v Va  - *>T  \ V..  , 


(2.8) 
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The  expressions  (2.3)  and  (2.9)  lead  to  a first  obser- 
T 

vation.  If  A - A is  small,  we  have  that  splitting  I 
results  in  a rather  fast  converging  method,  whereas  in 
splitting  II  the  convergence  behaviour  is  limited  by  a 
factor  A.  In  this  case  one  might  expect  splitting  I to  be 

>p  . i 

the  most  efficient  one.  The  assumption  A — A is  start*  > nor 

m 

unreasonable,  since  from  (?.l)  the  matrix  A - A arises 

T 

from  the  first,  order  term,  and  thus  we  have  A - A » 0(h) 

T 

compared  to  A + A . 

T 

One  might  however  put  tlve  question  what  happens  if  A - A 
is  not  very  small  ? 

From  matrix  theory  it  is  known  that  all  the  eigenvalues  X * 

T — 1 T " 

of  (A4A  ) (A-A  ) are  purely  imaginary  and  let  us 

define  as  the  maximum  absolute  value  of  these  eigen- 

\ 

values.  Then  for  splitting  I,  A ^ is  a measure  the  speed 

of  convergence,  while  for  splitting  II  the  convergence 
factor  behaves  like  ‘ | « Thus  we  conclude  that  for 

\ <1  ^ ^ , splitting  I leads  to  a divergent  process, 

while  splitting  II  still  yields  a converging  process. 

More  explicitly,  it  follows  that  for  - 1 c c VI s , splitting 
II  will  be  the  faster  one.  Therefore  the  question  arises 
whether  the  convergence  can  be  influenced  by  other  splittings 
of  the  matrix  A. 

Consider  now  the  following  class  of  splittings 

ft;  t>A  ( (W  P ^ (2.10) 

where  ^ i o is  an  arbitrary  real  constant.  The  splitting 
(2.10)  results  in  the  iterative  method 

, l - (2.11) 

And,  for  V*-  - -v  , we  have 

The  eigenvalues  of  the  matrix 

("  11  ‘ !(#. »■>)-’(» -«0 


-4- 


are  related  to  the  purely  imaginary  eigenvalues 
(A+AT)71(A-AT)  : 


For  wova  | A k | it  follows  that  the  choice 

' . l l X1 

-■  1 4 l 

leads  to  the  expression 

v ^ n 

5 >*-*«  A ¥ X'wa* 


V ; 1 - 1 X 

' J ia  2a 


A}  of 

(2.13) 

(2.14) 

(2.15) 


where  | ^ - | . 

Prom  (2.15)  it  follows  that  we  may  expect  convergence 
always,  if  o( . ot  t , since  ^ 1 . It  follows  imme- 

I 

diately  from  (i.  14)  that  If  A->A  , we  have  <>(„  t-  V » oince 
XK>A(>  0 * This  agrees  with  our  earlier  observation  that  the 

splitting  1 (2.2),  which  results  from  is  optimal  for 

almost  symmetric  matrices.  The  splitting  II  (2.3)  results 
from  (2.10)  if  we  choose  4 , and  this  choice  is  opti- 

:nal  if  * » 

The  next  observation  is  that  for  systems  that  arise  from 
problems  like  (2,1 ),  we  have  in  general  A r 0(V whore 
h is  a measure  for  the  gridsize  over  which  is  discretised. 
Prom  the  definition  of  i*  0ft  , it  follows  that  - 0(V*‘  ) 

This  indicates  that  the  proper  choice  of  c< 0^t  yields 
considerably  fiister  convergence  of  the  corresponding 
Iterative  method  (2.1 l).  The  last  observation  is  that  the 
eigenvalues  V-  are  all  situated  on  the  straight  line  X.  ~ 
in  symmetric  positions  to  the  X-axis.  One  could  use  this 
fact  in  order  to  follow  Manteuf fel ’ s ideas  for  acceleration 
of  the  Iterative  method  (2.11). 


‘i . Outer  aal  inner  iterations 


In  each  step  of  the  iterative  method  (2.11)  a linear 
system  of  the  form 

:!»'  (3.1) 

has  to  be  solved. 

Such  a step  will  be  called  an  outer  iteration  step.  If  the 
linear  system  arising  in  each  outer  iteration  step, 

is  solved  by  an  iterative  method,  then  the  steps  of  this 

d c « i 'bed  vjV  vi 

method  will  be  msarfcMtmfl  asVlnnor  iteration  steps. 

In  this  secti-n  we  will  consider  iterative  methods,  for 
solving  (3,l),  that  are  based  on  regular  splittings  /t> / 

fw  n 1 - ic  - n (5,?) 

For  these  regular  splittings^holda  U "5.0  and  R ° 
and  moreover,  if  R.  i 0 , 

0 c Ic  1 ^ + ' V ' (3.3) 

One  might  follow  different  strategies  in  the  outer-inner- 
iteration  process,  the  most  extremal  strategies  are  consi- 
dered here  in  some  detail. 

S t ra  t e gv  1 . The  equation  (?.l)  is  solved  at  each  outer— item— 
tlonstep  accurately,  that  is,  with  an  accuracy 
less  than  or  equal  to  the  desired  accuracy  in  the 
final  solution  of  Ax-b. 

strategy  ?,  Only  one  step  of  the  inneriteration  process  to 

solve  (3.1)  is  performed  at  each  outo r-itera tion 
step. 

It  might  be  expected  that  strategy  2 lends  to  an  Increase  in 
the  number  of  outer-iterations,  needed  to  reach  a certain 
accuracy,  as  compared  to  strategy  1.  From  a practical  point 
of  view  however,  It  Is  interesting  to  investigate  whether  the 
total  number  of  innerlterations  decreases. 


In  the  numerical  experiments  (section  •! ) it  whs  observed 
that  stretchy  was  the  cheapest  one.  It  was  even  observed 
that  strategy  ? needed  sometimes  1 e a a ou  to  r t te  rs  t ions  , which 
la  somehow  in  contrast  with  the  oonse *-va tionlaw  of  trouble. 
This  effect  will  be  more  or  leas  explained  here  for  the  case 
oi  * . For  a more  general  choice  of  ot  ■>  o , the  effects 

can  be  explained  similarly. 

For  strategy  2 it  follows  from  (i..’5  that 

Ic  \0tl  (5.4) 

and , i f we  set  " ' ' * ‘ ^ 

wo  havo 


R-'K-(A-fA)  ami 

w i,iM  [ nT*  u • co i f*1 ) J v’ t- 

or,  more  conveniently 


x - x 


then 


(?.s) 


V*.  • 


' (iWat\  \ u ’ (■  a nvy][\,t  (5.6) 


We  now  assume,  in  a very  rough  and  not  precise  way  that 


Vr  >•»  ( A » o1  ^ 1 

and  from  (5.0  we  have  r>  c.  I , 

Kquation  ( 5,o)  can  now  be  rewritten  in 

V-.. 


(5.7) 


' >M  *(0  . P’ 1 )‘(t>  0^ \\;  ( 5.0) 

For  the  eigenvalues  >5  of  the  matrix  l,  5\(#  d'tl  - il’ 5 

we  have  that 

1 iv"‘h'\  ( 5.9) 


c\ 


For  strategy  1,  we  have  iv,  1 , and  thus  A ■ . It  may  be 

expected  that  strategy  2 needs  less  outerl terations  if 

the  solution  of 

(5.10) 


bo  defined  ns 

r 


Let  therefore  \ 

■*o<< 

From  simple  cal cul a t Iona , it  follows  that 

(5.1D 

As  in  general  v,  will  be  very  near  to  1 , we  have  \2  I . This 

V\lt#  I 

oxpl/timi  fioniffhow  that  In  nt  ratify  ? Jx outer!  te  rat  iona  aro 


necessary  If  \. 
for  A < t 


K 

■r  ^ . , strategy  1 needs 

often  far  more  inner! tern t ions. 


I . However,  it.  should  be  stressed  ft ft—  , that 

i 

outeriterationa  indeed,  but 


Numerical  experiments. 


The  matrices  in  the  experiments  can  be  considered  as  to 
IT  t C > v\  i\  rnj 

jans wmmm  the  five-point  discretisation  of  the  partial 
dif ferential  equation 

-o  , (4.1) 

whore  U is  discretised  by  a central  difference  formula. 
The  equation  (4.1 ) is  given  over  a rectangular  region,  and 
u is  equal  to  a given  function  along  the  boundaries.  The 
matrix  arising  thus,  has  the  structure 


The  elements  of  the  diagonals  are  denoted  by  a^,  b^,  c^,  d^ 
and  e.,  where  i is  counted  rowwise.  As  a regular  splitting  for 
A+A  we  chooseJE  the  incomplete  decompoai tion  that  arises  in  the 
ICCG(3)-method  /4 /•  In  both  the  examples,  (4.1)  has  been  dis- 
cretised over  a rectangular  grid  with  10  meshpoints  in  the 
x-direction  as  well  as  in  the  y-direction.  ThiB  yielded  a 
matrix  A of  order  900  and  with  half  bandwith  30. 


Example  1 . In  the  first  example  we  demonstrate  some  of  the 

A/ *»C  < • kl-'  IV. 

effects  BWUftOOrtijtmlajt . Therefore  p , the  constant  in 

differential  equation  (4.1)  has  been  chosen  such,  that  the 

non-zero  values  of  the  elements  in  (4.2)  were: 

a^  =■  — 1 • i b^  = — 1.2S,  c^  4.  • = — «7b  * = —1. 

In  table  I the  results  for  different  strategies  are  listed. 

The  number  of  outeriterations  and  the  total  number  of  inner- 

iterations  to  reach  a certain  precision  are  represented. 

The  precision  was  estimated  by  Vj*»<  J X- . - X.  ( . | , where  X--  is 

the  i-th  element  of  the  outeriterationvoctor  x.  . From 

rather  rough  calculations  it  followed  that  /\^(the  max,  of 

T — 1 T 

the  absolute  values  of  the  eigenvalues  of  (A+A  ) (A-A  ) has 
a value  of  approximately  1.7. 
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Table  I.  Results  for  different  strategies 


Explanation  <a»  table  I: 

The  numbers  I up  to  VI  stand  for  the  different  strategies. 

The  number  of  outeriterations  is  counted  by  A and  B repre- 
sents the  total  number  of  inneriterationsteps.  The  strategies 


were : 


Is  (X;  ■(  'j  For  each  outeriterationstop,  the  inneriteration- 
process was  performed  in  an  accurate  way  (with  a 
maximum  of  10  steps').  This  is  strategy  1 (cor.v,  for 
A 

II:  ((in').  The  inneriterationprocess  was  stopped  a3  soon  as 
the  residual  was  less  than  1CK  of  the  initial  residual 
at  each  outeri terationstep  ( in  ||  ||  -norm). 

Ills  (j/-  ')  The  inneriterationprocess  was  stopped  a3  soon  a3 
the  residual  was  less  than  20 of  the  initial  residual. 
IV:  (ai  z -0  Strategy  2, 

V : \ } Strategy  2. 

VI:  (oi r l ) Strategy  1.  (convergence  for  Al_So<.<^  ) 

A few  additional  remarks: 

1.  As  X 3 Vi  it  could  be  expected  that  strategy  I resulted 
in  a convergent  process  and  strategy  VI  in  a divergent 
process . 

2.  It  is  surprising  that  strategy  V ( pi;  ^ ) leads  to  a 
convergent  process. 

3.  The  ideas  in  section  3 are  underlined  by  the  results  of 
strategies  I,  II,  III  and  IV. 

Example  2.  In  order  to  demonstrate  some  effects  for  a smaller 
^ X » (^  was  chosen  such  that  the  following  values  for  the 
nonzero  elements  in  (4.2)  resulted: 
a^=>— 1 . , b^»— 1.1  , c^=»4.  , d^=«— 0.9  , e^«—  1. 

From  rough  computation,  it  was  estimated  that  (cl). 

In  table  II  the  results  are  given  for  strategy  1 (accurate 
inneri tera tion)  and  strategy  2 (l-step  inneriteration) , both 
for  the  choice  in  (2,11). 


% Final  remarks 


It  was  observed  in  practical  situations  that  the  choices 
of  = or  o(~  1 in  the  iterative  method  (2.11),  in  combi- 

nation with  a 1-step  strategy  for  the  inneriteration  often 
lead  to  a fairly  efficient  and  easy  to  program  method. 

Only  the  inneriteration  based  on  an  incomplete  choleski- 
f actorization  (see  ICCG(})  in  /4/)  has  been  considered. 

No  attempts  have  been  4m*  to  accelerate  (2.11),  nor  has  kt 
b -en  tried  to  estimate  in  order  to  choose  an  optimum 

value  t»!  . It  should  be  mentioned  that  for  strategy  1, 

the  eigenvalues  ^ } of  the  iterationmatrix  are  also  situated 
on  a line  parallel  to  the  X-axis  in  the  complex  plane,  while 
in  strategy  2 the  eigenvalues  & • are  not  situated  on  such 

J 

a line,  but  in  the  neighbourhood  of  it. 
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Int  roduct ion 


Useful  variants  of  the  Lanczos  scheme  for  the  determination  of 

eigenvalues  of  large  symmetric  matrices  have  been  developed  in  the 

oast  few  years  (Paige  [6],  Golub  [7],  Lewis  [11]).  Symmetry  of  the 

matrices  is  essential  in  the  Lanczos  method.  Some  other  eigenvalue 

problems  can  be  reduced  to  symmetric  problems  after  some  preliminary 

T 

work  e.g.  if  the  matrix  A is  skew-symmetric  (A=-A  ) the  scheme  can  be 

2 

applied  to  the  matrix  A”,  which  involves  twice  as  much  computational 
work  (Cline  ; 1 2 I , Lewis  [111,  Platzman  [ 13]). 

Another  important  class  of  problems  is  concerned  with  the  determination 
of  eigenvalues  of  the  productmatrix  CB  where  C and  B are  symmetric 
matrices  and  one  of  them,  say  B,  is  positive  definite  as  well.  A common 

way  to  solve  this  problem  with  the  Lanczos  scheme  is  to  construct  first 

. . T T 

a Choleskt  decomposition  B=I.L  and  then  to  apply  the  scheme  to  L CL 

which  has  eigenvalues  identical  to  those  of  CB  (Golub  T7]).  Since  Lanczos- 

schemes  are  specially  attractive  for  sparse  matrices,  a disadvantage  of 

this  approach  might  be  a loss  of  sparsity  in  the  Choleski  decomposition. 

In  section  1 of  this  report  a generalized  Lanczos  scheme  is  proposed 
that  applies  directly  to  matrices  A whether  they  are  symmetric  or 
skew-symmetric,  and  to  produetmatrices  CB  where  C is  either  symmetric 
or  skew-symmetric  and  B is  symmetric  positive  definite. 

The  matrices  A,  B and  C do  not  have  to  be  represented  in  the  usual  way 
as  two-dimensional  arrays  of  numbers,  but  as  rules  to  compute  the 
products  Ax,  Bx,  and  Cx  for  any  given  x.  This  allows  us  to  take  full 
advantage  of  any  sparsity  structure. 

Lanczos  schemes  yield  approximations  for  the  eigenvalues  of  the  given 
eigenvalueproblem.  One  of  the  main  difficulties  is  how  to  distinguish 
between  good  and  bad  approximations,  since  both  are  generally  present 
(Paige  [61,  Parlett  and  Kahan  T 2 3 ) . In  section  2 an  algorithm  is  proposed 
to  determine  the  good  approximations  and  to  remove  the  bad  ones.  It 
should  be  mentioned  here  that  any  multiplicity  of  an  eigenvalue  of  the 
matrix  can  not  be  detected.  A multiplet,  if  there  is  one,  will  be 
represented  by  only  one  single  eigenvalue;  this  problem  is  peculiar  to 
Lanczos  schemes  (Kahan  and  Parlett  [2]). 


In  section  3 an  implementation  of  the  algorithm  of  section  2 is  given. 
Numerical  examples  for  the  algorithms  of  both  sections  I and  3 are 
presented  in  section  4. 

We  did  not  consider  in  detail  the  problem  of  determining  of  eigenvectors. 
In  section  5 we  summarize  the  main  results  of  Kalian  and  Parlett  T 2l  as 
well  as  some  of  our  numerical  results. 

Fortran  subroutines  for  both  the  generalized  Lanczos  scheme  and  the 
detection  of  good  eigenvalue  approximations  are  covered  in  the  final 
section  of  this  report. 
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A generalized  Lanezos  scheme 

In  this  section  we  describe  Lanezos  schemes  that  apply  to  skew- 
symmetric  matrices  or  product-matrices  CB,  where  B is  symmetric 
positive  definite  and  C is  either  symmetric  or  skew-symmetric 
(all  the  matrices  A,  B and  C will  be  of  order  n) . It  should  be 
mentioned  here  that  the  eigenvalues  of  CB  are  identical  to  those 
of  BC. 

Eigenvalue  problems  Cx=ABx  can  be  reduced  to  either  of  these 
forms . 

The  algorithms  are  closely  related  to  an  algorithm  published  by 
Widlund  [51  for  the  solution  of  non-symmetric  linear  systems. 

The  eigenvalue  problems  of  a skew— symmetric  matrix  A can  be 
reduced  to  the  eigenvalue  problem  of  a symmetric  matrix  by 
squaring  the  matrix:  A".  This  is  not  necessary  in  our  formulation, 

V 


1 


Pet  nil  turn  o t tln>  gene  i a 1 i :*ed  l.auczos  scheme 


I.«*C  A bo  of  tin-  form  A-CR,  where  11  is  symmetric  positive  definite 
and  C is  either  symmetric  or  skew-symmet  r ic . 

Choose  an  arbitrary  veer  or  v^,  with  iv  j »v  | • •1IK'  form  u “Av  . 

Kows  {v.},  lot.),  l|l.l  and  (y.l  are  then  generated  by 


“j  • <vj  • Avi'» 


W . “ U . - Cl  , V . 

j .1  J .1 


Vi  “ (”j  • ”i,» 


"j.l  - ''jM 


V . • - — \V  , 

1+1  \jH  1 


1 i “ Av  . , - (1 . v . 
, t- 1 _t+  1 j+  I j 


. ,m 


(as  fas  as  > . •fO, 
see  not  e 1 1 ' 


where  (x,y)  -(x,Wy)  , with  11  symmetric  and  positive  definite, 
and  t - I it  C“C * 


t — I if  O-C 


V 


(see  also  note  .’) 


I he  constants  on,  (l.  and  y.  define  a tridiagonal  matrix  T : 

ill  m 


T - 

in 


ft  j tV, 

Y , it , 


> it 

m m 


Not  i 


It  in  some  stage  >j“0,  then  one  ran  either  restart  with  a new 
v,  or  proceed  with  \j  replaced  by  some  small  constant.  In 
practice  the  situation  Yj-o  occurs  very  incident al Iv.  In  our 
implementation  such  a is  replaced  bv  a small  constant. 


Not«'  : For  It- 1 and  t-l  we  have  the  original  l.anozos  scheme  as 
defined  hv  Paige. 


Theorem 


T T 

We  assume  that  either  OC  (x-1)  or  O- C (x— 1)  holds  and  that 

B is  a positive  definite  symmetric  matrix  and  A-CB,  then  the 

generalized  Lanczos  scheme  applied  to  A generates  a tridiagonal 

matrix  T , where  limit-values  of  the  eigenvalue  of  T for 
m m 

increasing  m,  should  be  equal  to  some  of  the  eigenvalues  of  A, 
but  they  may  differ  by  a certain  amount  depending  on  the  precision 
of  computation. 


Proof 

T 

i)  For  OC  and  B=l,  the  result  is  well  known  (Paige  T6],  Golub  [7]). 

T 

ii)  For  C=-C  and  B=I  the  proof  is  as  follows: 

It  is  only  necessary  to  establish  that  the  generated  row 

lv,  },  . is  an  orthonormal  row.  The  nroof  is  by  induction, 

k k=  1 , . . m • 

Let  {v.  , .be  an  orthonormal  row. 

k k= 1 , , . . , j 

Then  we  have,  for  v.  , the  relation: 

J + > 


Y-  ,v . . 1 *=  Cv,  - B.v.  , - a.v,  , 

'j  + l J+l  j j j-l  J j ’ 

where  we  assume  that  y.  ,+0,  since  in  that  case  the  recurrence 

1 j + l 1 

relation  terminates. 

For  k<j-l : 

(YJ+ivj+rV  - (Cvj  " Pjvj^ 


- (v.,cvk) 


" (vi’Yk+ivk+i  + pkVi  + VkJ 


» Q 
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For  k= i- 1 : 


(Yj*ivi*,,vi-i>  ' (cY,i-i>  - 


“YYl’  - B. 


Since  3.  - -y.  « - (y.v.,v.)  =■  - (Cv.  ,,v.)  =»  (Cv.,v.  .) 
J J J J J J-l  J J 


it  follows  that  (y.  ,v.  ,,v.  ,)  = 0 

'j  + l j + i’  j-l 


For  k= i : 


‘W'-V  " “YV  ' “j  ’ 0 


Finally,  we  have 


= -~o — (Av.  - 0.v.  - a.v. ,Av.  - 3.v.  - a.v.) 

J + > J+>  v2  J J J-l  J J J J J-l  J J 

Tj+1 


V - - (u.  - a.v.,u.  - a.v.) 
2 J J J J J J 


— (wj.w j) 

Yj  + 1 


Thus  the  row  {v.  . ...  is  an  orthonorraal  row. 

k k=  I , . . . , j + 1 


- ; - 


iii)  When  OC  and  B is  symmetric  positive  definite,  B can  be 
T 

written  as  B=LL  , where  L is  lowertriangular . 

• T 

Since  the  eigenvalues  of  CB  are  equal  to  those  of  L CL,  the 

• • i T 

original  Lanczos  scheme  might  be  applied  to  L CL  (with  normal 

inner-product  ( , ) ). 

In  this  case  we  then  have  the  special  relation 


Oj  •>  (v^,L  CLVj ) 


uj+]  = (L  CLvj+1  - Pj+Jv.) 


it  follows  that 


Lu.  = LL  CLv,  , - P.  Lv . 

j+1  j+I  Kj+I  j 


If  we  replace  x by  L x,  this  equation  can  be  rewritten: 


, T„. . T, 


LL  uj+]  - LL  CLL  vj+]  - Pj+1LL  v. 


u . = CBv . - p . v . 

j+1  j+1  j + 1 j 


A5j.i  - “j.i’j 


The  other  Lanczos  relations  follow  from 


cx . = (L  CLv . , v . ) 

J J J 

= (LTCLLTv.,LTv.)  - (CBVj , BVj ) 
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(W..Vj) 


T T 

(L  w.,L  w.) 


- (Bw . ,w . ) 
J J 


(W.,w.) 


B 


The  relations 

w.  “ u.  - cx.v. 
J J J J 


and 


1 


are  obvious. 

The  vectors  w.,  v.  and  u.  produce  the  desired  result. 
J J J 


iv)  The  remaining  case  A=CB,  where  O- Cl  and  B is  symmetric 
positive  definite,  follows  from  the  previous  ones  (with 
T-l).  //. 


Remarks 


i . 

2 


T 

If  O-C  , we  have  ct.^O  for  all  j. 

J 

The  above  theorem  allows  the  computation  of  the  eigenvalues  of 

CB,  which  are  equal  to  those  of  BC,  without  the  explicit  need 
T 

for  an  LL  -factorization  of  the  matrix  B.  This  makes  the  new 

schemes  very  attractive,  especially  if  B has  a sparse  structure. 

However,  it  should  be  noted  that  eigenvectors  cannot  be  computed 

. T 

by  these  schemes  directly,  since  then  an  LL  -factorization  is 
required  for  a proper  transformation. 

Eigenvectors  may  be  computed  by  a Raleigh-quotient  iteration 
scheme  C 1 1 , once  one  has  a (fast)  solver  for  systems  like  Bx“y. 

For  sparse  matrices  B,  for  which  fast  direct  or  iterative 
solution  schemes  exist,  this  Lanczos  scheme  can  also  be  used  for 
determining  eigenvalues  of  Cx-\Bx,  via  B *Cx«\x.  The  scheme  should 
be  applied  to  CB  ' which  has  identical  eigenvalues. 
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3.  We  should  like  to  mention  briefly  certain  aspects  of  programming. 

For  the  generalized  problem  the  adapted  schemes  O.l)  require 

only  one  extra  matrix-vector  multiplication  and  only  one 

additional  vector  to  store  Bw..  Remember  that  Bv.  can  be 

J J 

computed  from 


4.  If  C is  skew-symmetric  (x=-l)  then  the  generated  matrices  T 

m 

are  also  skew-symmetric.  Eigenvalues  of  a tridiagonal  skew- 
symmetric  matrix  can  be  computed  as  follows; 

the  matrix  iT  is  Hermitian  and  has  real  eigenvalues.  Since  in 
in 

the  computation  of  the  eigenvalues  with  Sturm-sequences , only 

2 

squares  of  off-diagonal  elements,  , are  involved,  these 
eigenvalues  can  be  computed  without  any  complex  computation. 
Once  the  eigenvalues  of  |T  |: 


b2  0 63 


,TJ  5 


have  been  computed,  they  should  be  multiplied  by  i so  that  they 

represent  the  eigenvalues  of  T . 

in 
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Monitoring  of  the  Lanczos  process 

In  the  single  vector  Lanczos  processes,  as  described  in  section  1, 

rows  of  mutually  orthonormal  vectors  are  generated. 

The  coefficients  a.,  8.  constitute  a tridiagonal  matrix  T the 
J J m 

eigenvalues  of  which  bear  some  relation  to  those  of  A. 

We  only  consider  here  the  symmetric  case  (x=l),  the  other  one 

(x"-l)  is  obvious. 


ci , 8i 


B2  a2  83 

83  a3 


8 , a . 8m 

m- 1 m- 1 m 1 

8 a 

m m 


For  B=I,  T is  related  to  A by 
m ' 


AV  = V T + 8 ^.v  x1  e* 
m mm  m+1  m+1  m 


where  V is  an  orthonormal  n*m  matrix,  in  which  the  v are  the 

m T 1 

columns,  and  e »(0,0, . . .0,1) , the  m-th  unitvector  (n  is  the  order 
m 

of  the  matrix  A) . 

The  relation  of  the  eigenvalue-problem  of  T^  to  the  eigenvalue- 
problem  of  A is  discussed  and  demonstrated  by  extensive  numerical 
experiments  reported  by  van  Kats  and  van  der  Vorst  [3], 
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2a.  The  eigenvalue  problem  T x=Xx. 

n; 


The  eigenvalues  of  T are  the  roots  of 
m 


(2.3) 


det (T  - A)  - 0 
m 


If  the  leading  k-th  order  principal  minor  of  T^,  is  denoted  by 
T , then  the  following  recurrence  relation  holds 


(2.4) 


det(Tk  - A)  - (c^  - A)det(Tk_!  - A)  - &£det(Tk_,  - A) 


If  we  define  det(Tg-A)=l  and  since  we  have  det (Tj-A)=Cij-Af  it  follows 
that  the  above  relationship  is  exactly  that  of  orthogonal  polynomials 


(2.5) 


pk(x)  = (^  - x)  pk_,(x)  - \Pk_2(x> 


(2.6) 


It  follows  that  the  zeros  of  p separate  the  zeros  of  p,  , (x) 

k k-1 

strict  sense 

(Wilkinson  ( 1 1)  . 


and  pk+j(x)  in  a strict  sense  if  none  of  the  3^  equals  zero 


By  analogy  then,  if  the  ordened  eigenvalues  of  T are  denoted  by 

. (k ) , R 

A . , we  have 

l * 

»‘k;»  < »<k>  < 

l-l  i l 


2b.  Recognition  of  limitvalues 

From  relation  (2.6)  it  follows,  that  the  extreme  eigenvalues  if 

T , for  increasing  m,  converge  strictly  monotonically . Since 

according  to  Paige  limitvalues  of  the  row  T are  equal  to  eigenvalues 

m 

of  the  original  matrix  A,  except  for  some  amounts  that  depend  on 
the  precision  of  computation,  this  property  may  be  used  for  an 
automatic  determination  of  the  extreme  eigenvalues.  However,  it  is 
evident  that  limit  sequences  can  also  be  recognized  for  internal 


25 


12 


eigenvalues,  since  the  strict  separation  relationship  (2.6) 
should  hold. 

(k) 

As  soon  as  relation  (2.6)  is  violated,  either  because  X)  equals 
one  of  the  or  xj^  or  is  outside  the  interval  [xf^j'\xj^  '^] 

we  know  that  at  least  in  the  precision  in  which  we  are  working  it 
is  not  possible  to  distinguish  between  X.  and  the  respective 
eigenvalue  of  T^_j.  Consequently  we  have  a limitvalue  and  thus  an 
approximate  eigenvalue  of  A.  Since  we  are  working  in  a finite 
precision  and  the  extreme  eigenvalues  are  bounded  by  the  extreme 
eigenvalues  of  A,  the  separation  relationship  (2,6)  will  be  violated 
sooner  or  later. 

In  practice  this  provides  us  with  an  excellent  means  of  recognizing 
limitvalues  automatically. 


As  soon  as  one  case  of  violation  has  been  encountered  one  measures 

how  much  the  respective  values,  say  b and  c,  differ  relatively. 

A value  e.  is  defined  as  follows 
be 


(2.7) 


abs (b-c) 
Cbc  l+abs(b) 


The  maximum  over  all  violations  will  be  taken  as  e.  This  e yields  an 
impression  of  the  relative  accuracy  with  which  all  the  eigenvalues 
T^  and  T^_ ^ have  been  computed. 

(N.B.:  this  is  not  to  be  confused  with  the  accuracy  of  the  Lanczos- 
process  itself). 

As  a rule  of  thumb  this  e is  multiplied  by  n (the  order  of  the 
original  matrix  A)  and  all  eigenvalues  of  and  T^  ^ which  differ 
by  less  than  n*£  will  be  taken  as  possible  limitvalues. 

In  the  next  section  this  will  be  stated  more  precisely. 
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KVSCAN;  implementation  of  the  monitoring  process 

The  monitoring  process  is  essentially  based  on  the  separation 
relationship,  as  described  in  the  previous  section.  This  requires 
the  computation  of  all  the  eigenvalues  of  two  succeeding  tndiagonal 

matrices  and  T^. 

At  the  first  stage  of  the  process  one  checks  to  see  whether  the 
separation  relationship  has  heen  violated.  It  is  well  known  that  in 
the  Lanczos-process  multiplots  of  eigenvalues  are  introduced  as  soon 
as  orthogonality  has  been  lost  P3. 

The  next  step  is  to  recognize  these  multiplets.  Each  multiset  will 
be  represented  by  one  single  eigenvalue  interval.  If  the  original 
matrix  A has  a multiplet  eventually,  this  will  not  be  recognized. 

After  the  eigenvalues  of  and  Tfc  have  been  scanned  for  multiplets, 

die  resulting  multiplet-f ree  rows  are  compared  in  order  to  determine 
those  limitvalues,  which  represent  eigenvalues  of  A. 

The  monitoring  process  will  be  described  in  detail  below. 


Step  l 


'Check  whether  the  eigenvalues  if0  of  \ separate  the  eigenvalues 
x(k-l)  o(  T in  a strict  sense.  If  some  X i is  outside 
the  interval  [X^"’ } ,x[k_’ ) 3 then  this  yields  a lowerbound  e’ 
for  the  highest  attainable  relative  precision  in  all  the 
eigenvalues,  and  we  define  e”=max  e\  where  the  maximum  is  taken 
over  all  violations.  If  no  violation  has  been  encountered  then 
e"  is  taken  to  be  2_£,  where  t is  the  number  of  digits  in 
floating  point  arithmetic. 

An  upperbound  for  the  relative  working  precision,  to  be  used 
in  the  following  steps,  is  estimated  by 


where  n is  the  order  of  the  matrix  A. 


frump 
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Step  2: 


(k_n  (k\ 

With  the  e resulting  from  step  I,  the  row  (X)  and  (X)  '} 

are  both  scanned  separately  for  multiplets.  As  soon  as  a 

multiple  value  has  been  discovered,  i.e.  two  values  are 

encountered  which  differ  relatively  by  less  than  e,  the 

eigenvalue  concerned  is  represented  as  an  interval  with  the 

smallest  value  of  the  multiplet  as  the  lowerbound  of  the 

interval  and  the  largest  one  of  the  upperbound.  If  successive 

eigenvalues  are  recognized  as  belonging  to  the  same  multiplet, 

this  may  lead  to  a larger  value  for  the  relative  precision 
c ( a^s  Cupperbound— lowerbound)  . 

l+abs( lowerbound)  ‘ 

Step  2 is  repeated  with  the  most  recent  value  of  e as  long  as 
£ increases. 


Step  3: 

From  step  2 two  rows  of  intervals  result,  representing 
eigenvalues  of  T^_ ^ and  respectively.  These  rows  are,  as 
far  as  possible,  multiplet-free  with  respect  to  c • 

For  each  interval  in  one  row  one  checks  to  see  whether  there 
is  an  interval  in  the  other  row  that  intersects  with  the  first 
one  or  is  at  a distance  of  less  than  e relatively.  If  one  of 
these  conditions  has  been  met,  a new  interval  is  chosen  as  the 
span  of  both.  The  length  of  the  new  interval  yields  a new 
value  for  e. 

Step  3 is  repeated  with  the  most  recent  value  of  e as  long  as  e 
increases . 

If  6 successive  intervals  in  this  process  are  encountered 
belonging  to  T^_ ^ or  T^  for  which  the  above  condition  does  not 
hold,  then  a hole  in  the  spectrum  is  assumed.  The  value  6 has 
been  chosen  from  numerical  experience  and  could  be  replaced  by 
any  better  value. 

Continuing  in  this  fashion,  step  3 delivers  one  single  row  of 
intervals,  each  of  which  may  be  considered  to  contain  a limitvalue. 

These  limitvalues  differ  only  from  the  eigenvalues  of  the  original 
matrix  A with  regard  to  the  degree  of  precision  in  which  we  are  computing. 
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In  some  situations  it  may  occur  that  step  3 yields  an  interval 
which  contains  no  eigenvalue  of  A.  However  in  such  cases  there  is 
an  eigenvalue  in  the  neighbourhood  of  the  interval.  In  these 
situations  it  is  common  for  the  process  to  yield  also  the  interval 
in  which  the  respective  eigenvalue  is  situated;  thus  two  very 
close  intervals  are  obtained. 

In  order  to  identify  both  intervals  as  representing  the  same 
eigenvalue,  it  is  advisably  to  apply  only  step  2 to  the  final  row 
with  a slightly  larger  value  for  e (I0*e,say). 


PStOG=,a‘,  LNCZOS  (OUTPUT  ) 

OI.iENS  ION  '30(2l),31(ci),y(21),aLDHfl(16)«B£TA(lfe) 
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Numerical  Experiments 


All  the  numerical  experiments  have  been  carried  out  on  the  CDC 
Cyber  73-28  of  ACCU.  The  relative  working  precision  is  48  bits 
(about  14  decimal  digits). 


4.  1 The  Bar-problem 

The  bar-problem  is  known  to  cause  difficulties  in  the  determination 
of  the  eigenvalues  at  the  lower  end  of  the  spectrum  T4,  3) 

The  40th  order  matrix  is  given  by: 

/ 

5 -4  1 

-4  6-4  1 

1-4  6-4  1 


0 1-46-4. 

1-4  6-4 

1 -4  3 

' ) 

With  EVSCAN  the  eigenvalues  are  determined  in  each  10th  step  of 
the  Lanczos-process  i.e.  T,.  , is  compared  to  T , . . Tn  table  1 

we  summarize  the  number  of  different  eigenvalue  intervals,  as 
delivered  by  EVSCAN  for  k=3  up  to  40.  Ve  note  that  when  k is  29, 

30  or  37  more  than  '40  eigenvalues  are  found.  If  we  follow  the 
advice  and  make  an  extra  scan,  as  described  in  section  2,  then  in 
general  the  number  of  eigenvalue  intervals  is  not  affected,  except 
when  k is  29,  30  or  37.  In  the  latter  cases  the  final  number  of 
eigenvalues  was  corrected  to  40  at  the  cost  of  some  slightly  larger 
intervals.  This  implies  that  there  have  been  migrating  eigenvalues 
very  close  to  an  eigenvalue. 
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From  the  experiments  it  follows  that  all  eigenvalues,  including 
those  at  the  lower  end  of  the  spectrum,  are  detected  by  F.VSCAN  at 
some  stage.  In  order  to  determine  the  eigenvalues  at  the  lower 
end  one  has  to  perform  a number  of  Lanczos  steps  larger  than  the 
order  of  the  matrix  (70  steps;  order  is  40). 

From  table  I it  can  be  seen  that  the  estimated  e,  scaled  with 
respect  to  the  machine-accuracy,  increases  slowly.  Since  it  is 
natural  to  relate  the  accuracy  to  the  order  of  the  respective 
tridiagonal  matrices  it  can  be  seen  from  the  fourth  column  of 
table  1 (t'/k)  that  the  eigenvalue  determination  with  this  process 
is  rather  stable. 


W i 1 k i nson ' s 


W2|  and  W2I 


Wilkinson  I I 1 has  introduced  two  classes  of  tridiagonal  matrices 

W*'  and  W„  . . 

2n+  1 2n+  I 


W_  , is  defined  bv  the  relations 
2n  H J 


ow  c n + 1 - i (i“l,...n+l),  tt.  *>  i - n - 1 ( i*=n+2 , . . . 2n+ 1 ) 

P.  = 1 (i  = 2 .... 2n+  1 ) 

and  W , by  the  relations 
2n+  1 

n.  » n + I - i (i=!,...n+l)  , ot.  = n + 1 - i ( i=n+2 , . . . 2n+ 1 ) 

P{  “ I (i-2,...2n+l) 


, T — 

The  matrices  and  W.,  have  been  used  for  our  tests. 


+• 

w, 


71 


10  I 
1 9 1 


8 I 


1 0 1 


1 8 1 


V. 


1 


9 


10 


10  1 


I 8 I 


1 0 1 


1 -10 


Table  11  g i vos  the  eigenvalues  of  W.,j  and  table  Til  those  of  W,j 

. ♦ 

As  (‘.in  be  seen  Iron)  table  11  some  of  the  eigenvalues  of  WP|  are 
Unite  ('lose.  With  the  precision  in  which  wo  are  eompnt  ing  we  cannot 
expect  to  detect  .’1  separate  eigenvalues;  at  least  A^  and  A ^ 
might  behave  as  multiple  eigenvalues. 

Since  the  error  in  the  l.anc/os-procoss  increases  slowly  (set1  4.0, 

*■  . ... 

W.,j  has  been  chosen  in  order  to  determine  which  eigenvalues  are 

recognized  as  being  distinct  in  several  stages  of  the  process. 

Also  we  may  get  some  impression  of  the  behaviour  of  the  l.anczos— 

process  for  (almost)  multiple  eigenvalues.  The  motivation  tor  the 

choice  of  W , ^ is  somewhat  different.  The  eigenvalues  of  W , ^ are 

equal  and  opposite  in  pairs.  It  is  well  known  that  the  power  method 

gives  slow  convergence  in  such  cases  I 1 ’] . Since  there  is  some 

relationship  between  the  powermethod  and  the  bane zos -method  it  might 

be  interesting  to  apply  the  I.anozos-proeess  to  this  matrix. 

With  respect  t o the  numerical  experiments  for  W,^,  the  following 
observations  have  been  made. 

- From  a selection  of  the  results,  as  represented  in  table  IVa,  it 
follows  that  in  no  case  all  .'1  eigenvalues  are  detected.  Only  for 
m .’0  eigenvalue  intervals  have  been  determined,  for  higher 
values  ot  m even  A^  and  A are  represented  by  one  eigenvalue 
interval.  This  explains  also  the  increase  in  t . 
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- A more  serious  defect  is  the  following.  With  almost  multiple 

eigenvalues,  which  are  distinct  within  our  accuracy  of  computation 

( X . „ , X,„  and  X,_,  X,,)  eigenvalue  intervals  are  delivered  which  are 
IV  I o 1 / lo 

in  between  both  eigenvalues.  These  eigenvalue  intervals  differ 
significantly  from  both  true  eigenvalues,  see  table  Va  (eigenvalue 

.8038 E+0  I and  .9210 E+0 1 ) and  compare  these  values  with 

those  in  table  II. 

- Table  Va-d  list  detailed  results  for  W0 ^ . 

With  respect  to  W f ^ it  is  observed  that  there  are  no  problems  in 
determining  the  eigenvalues,  except  with  respect  to  the  speed  of 
convergence.  For  m=16  (see  table  IVb)  no  eigenvalues  have  been 
detected  automatically. 

Detailed  results  are  listed  in  table  Vla-c. 


4 . 3 Pathologically  clustered  eigenvalues 


It  is  well-known  that  the  convergence  properties  depend  highly  on 
the  relative  clustering  of  the  eigenvalues.  Therefore  we  have 
constructed  a matrix  with  a cluster  of  eigenvalues  at  each  end  of 
the  spectrum  and  one  single  eigenvalue  in  between  both  clusters. 

The  matrix  A chosen  was  a 40th  order  diagonal  matrix  with  diagonal 

elements:  1.001,  1.002,  1.019,  2.000, 

3.021,  3.022,  3.039,  3.040. 

After  10,  20,  30  steps  of  the  I.anczos-process  only  the  eigenvalue 
2.0  is  determined  automatically  (no  convergence  at  the  lower  and 
upper  end  of  the  spectrum). 

After  40  Lanczos  steps  convergence  was  signalled  at  both  lower  and 
upper  end  of  the  spectrum  but  the  second  eigenvalue  of  A (1.002) 
was  not  found.  For  results  see  table  Ilia. 

After  45  Lanczos  steps  all  eigenvalues  have  been  determined,  see 


The  Lanczos-method  is  proposed  usually  for  the  determination  of 
the  extreme  eigenvalues  of  sparse  matrices.  We  have  used  the  Paige 
style  Lanczos  algorithm  to  compute  some  of  the  extreme  eigenvalues 
of  a symmetric  full  matrix  of  high  order  (n=*519). 

The  matrix,  used  in  this  example,  originates  from  a nuclear  shell- 
model  calculation  18  I. 

In  such  a calculation  one  computes  the  matrix  elements  of  a given 
one  plus  two-body  interaction  Hamiltonian  in  a set  of  j-j  coupled 
many  particle  basis-states.  After  diagonal isat ion  one  obtains  the 
energies  and  the  wave-functions  of  the  systems.  In  the  present  case 
the  basis  chosen  is  approximate  for  the  description  of  nuclear 
states  in  with  zero  angular  momentum  and  positive  parity. 

The  order  of  the  matrix  is  510;  it  contains  about  58%  zero-valued 
elements.  No  advantage  has  been  taken  of  the  zero  values  which  are 
distributed  in  rather  an  irregular  way. 

This  matrix  has  well  separated  eigenvalues  (no  multiplets),  which 
are  distributed  over  the  interval  (-160.0,  -180.5). 

In  general  the  eigenvalue  problem  for  full  symmetric  matrices  is 
solved  by  the  Householder  method  [ll. 

For  matrices  which  cannot  be  stored  in  fast  core,  this  process  is 

complicated  to  programme.  Since  the  CP-time  used  is  roughly 

3 

proportional  to  1/3  n (n  is  the  order  of  the  matrix),  the  Lanczos- 
process  for  the  determination  of  the  extreme  eigenvalues  is 
advantageous  if  less  than  1/3  n iterations  are  required. 

Other  practical  advantages  of  the  Lanczos-process  in  this  case 
are  that  it  is  easy  to  restart  and  easy  to  programme . 

In  table  Vllta-b  the  results  are  listed  for  40  and  60  Lanczos  steps 
respectively;  in  the  latter  case  the  scanning  process  has  been 
performed  with  a larger  r too  (table  VI  lie). 

A larger  f has  been  chosen  since  the  process  for  a large  Cull  matrix 
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is  generally  expensive  and  one  wants  to  extract  as  much  information 
as  possible  from  the  iteration-steps  performed.  The  possibility 
of  scanning  with  a larger  e has  not  been  included  in  EVSCAN  but 
could  be  with  a minor  modification. 

For  results  see  table  Vllla-c.  For  nv*40,  convergence  at  the  lower 
end  of  the  spectrum  is  detected,  for  m"60  convergence  is  detected 
at  both  the  lower  and  the  upper  end. 


4 . 5 A large  sparse  matrix 


For  the  matrix  used  here,  we  have  chosen  the  modified  Laplace 
problem  as  describe  in  [141. 

The  matrix  A results  from  five  point  discretisation  of  Au=0  over 

the  square  region  0<x,  yN 1 . The  boundary  conditions  are  ~ = 0 
— 3u 

for  x=0  x=l  and  y=l,  and  u=I  for  y=0. 

This  equation  was  diseretised  over  a rectangular  grid  with 
meshspacings  h = and  h = ~ , thus  yielding  a matrix  of  order 
992.  For  this  matrix  the  ICC0(3)  decomposition 


A - L^]  + R3 

is  constructed  (see  [9]).  The  eigenvalues  of  L3  AL3  have  to  be 
computed.  They  are  very  strongly  clustered  around  the  value  1.0. 
Also  at  the  upperside  of  the  spectrum  the  eigenvalue  distribution 
is  very  dense.  The  largest  eigenvalue  is  1.17  and  the  smallest  one 
is  close  to  0.0. 


In  table  IX  the  results  (number  of  eigenvalues  and  scaled  e)  are 
listed  for  several  stages  of  the  I.anczos-process  (m  denotes  the 
number  of  steps). 

In  figure  I the  following  quantities  are  represented  as  a function  of 

- the  total  number  of  eigenvalues,  detected  by  EVSCAN  (x) 

- the  number  of  eigenvalues  at  the  lower  end  of  the  spectrum  (0). 

This  number  was  fairly  well  represented  by  the  parameter  NT1IV 
(see  the  description  of  EVSCAN). 

- the  number  of  eigenvalues  at  the  upper  end  of  the  spectrum  (+) . 


m 
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4.6  A generalized  eigenvalue  problem  CBx=*Xx 

In  order  to  demonstrate  the  applicability  of  the  generalized 

Lanczos-scheme , as  described  in  section  1,  the  following  problem 

has  been  chosen.  Some  of  the  eigenvalues  of  CB  will  be  computed, 

where  B is  the  5 point  finite  difference  approximation  of  the 

Poisson-operator  A over  a square  10*10  grid  and  C is  the  central 

g 

difference  approximation  of  the  operator  ^ on  the  same  grid. 

The  matrices  look  like  this 


C = 


0 +0.  I 

\ \ 

-0.1  x v 
X X 

\ \ \ 

X \ * 

\ x+0. 1 

-0.  1 0 


T0  +0.1 

i-^\ 

! XS  s\ 

\ \ * 
x \0-' 
0.  1 0 


blocksize  is  10, 

10  blocks  along  the 
diagonal : order  of 
C is  100 


B = 


4 -1 


x 


X 

\ \ X 

X X 


-I 


x -1 
-1  ' 4 


4 -1 

X 'x 

“l  \ ' 

X 

X X 


x -1 

-1  S 4 


same  structure 
as  C. 


The  results  are  listed  in  table  X for  m“IS,  30,  60,  where  m is  the 
order  of  the  tridiagonal  matrix,  generated  in  the  generalized 
Lanczos-scheme . 
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5 . Notes  on  eigenvectors 


5 . I Theoretical  Aspects 


Before  we  mention  some  results  of  our  numerical  experiments, 
relevant  results  of  Kahan  and  Parlett  [2]are  summarized. 

The  results  of  the  Paige  style  Lanczos  scheme  for  a symmetric  or 
skew-symmetric  n-th  order  matrix  A can  be  written  as: 


AV  = V T + 6 v e 
k k k Pk+1  k+1  k 


where  e =(0,0,0, ... ,0, 1)  , the  k-th  unitvector. 


V is  formed  by  the  columns  v. 
k l 


1 3 1 1( 

1 I I • t • y IV  « 


For  any  y and  normalized  vector  x,  the  quantity  | |Ax-yx| | bounds 
the  error  between  y and  some  eigenvalue  X of  A,  If  y is  an 
eigenvalue  of  T^  and  y the  associated  eigenvector  then 


| X - y | < | | Ax  - yx | 


lAV.y  - yV.y| 


•ky  " VkV' 


lyJ 


Thus  the  error  in  this  computed  eigenvalue  is  bounded  by  |Pjc+jl 
times  |y^|  and  if  y^,  which  is  the  last  component  of  the  vector  y, 
is  small,  then  the. bound  may  be  sufficiently  small  even  though 
8k+|  is  of  moderate  size.  From  this  analysis  it  follows  that  the 
more  accurate  approximations  to  eigenvalues  may  be  those  eigenvalues 
of  T^  whose  associated  eigenvectors  have  rapidly  dwindling  components. 


I 
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5.2  Numerical  results 


Eigenvectors  of  the  original  matrix  A have  been  computed  in  the 
following  way. 

With  EVSCAN  eigenvalue  intervals  of  the  final  tridiagonal  matrix 

have  been  determined.  Since  TST11RM  (Eispack  C 1 0 1 > requires 

intervals,  these  intervals  could  be  supplied  directly.  Since  an 

eigenvalue  of  T could  be  equal  to  an  upperbound  or  a lowerbound 

of  an  interval,  which  is  not  permitted  by  TSTURM,  the  intervals 

have  been  made  slightly  larger.  The  eigenvectors  of  T have  to  be 

back transformed  by  V to  eigenvectors  of  A. 

k T 

For  productmat rices  an  extra  LL  -decomposition  of  the  matrix  B 
(see  section  1)  is  required  for  backtransformation,  thus  nullifying 
the  advantages  of  the  generalized  Lanczos— scheme. 

In  the  experiments  we  demonstrate  the  behaviour  of  the  last 

components  of  a normed  eigenvector  of  T . To  this  end,  the  matrix 

+ . k 
W , | (see  section  4)  has  been  chosen. 


5.2.1  For  k=!3  no  eigenvalues  could  be  detected  automatically  by  EVSCAN; 
this  is  reflected  by  the  behaviour  of  the  last  components  of  some 
eigenvectors . 

A Lanczos-approximation  for  the  eigenvalue  2.130209219163  was: 

2.  13327 The  last  4 components  of  the  corresponding 

eigenvector  of  T are: 

0.31  , -0.35  , -0.013  , 0.31  , 

These  components  also  indicate  that  no  convergence  had  occurred. 

Also  for  k-13,  the  largest  eigenvalue  of  W*  , 10.746194182903  has 
been  approximated  by  10.746194182902.  The  better  convergence  is 
reflected  by  the  last  4 components  of  the  corresponding  eigenvector 
of  T : 


-.0041  , -.0010  , -.00013  , -.000013 


evident  that  limit  sequences  can 


5.2.2  For  k=16,  EVSCAN  detected  a number  of  eigenvalue  intervals 
(see  section  4) . 

A Lanczos-approximation  X for  the  eigenvalue  10.746194182903  was 
given  by  the  same  value.  The  last  4 components  of  the  corresponding 
eigenvector  of  are: 

. 1 3E-04  , . 27E-06  , .57E-09  and  .59E-10. 

If  we  denote  x as  the  backtransf ormed  eigenvector  of  W2],  Chen  we 
t |W^  x -Ax  || 

have  = . 35E- 10 


In  this  case  EVSCAN  yields  an  eigenvalue  interval  which  does  not 

contain  an  eigenvalue  of  W,,  (though  there  is  one  close  by). 

The  eigenvalue  7.003952209528  was  approximated  by:  7.003952209098. 

In  this  case  the  last  4 components  of  the  eigenvector  of  T , are: 

1 o 

. 21E-0I  , . 89E-03  , .48E-05  and  .15E-05. 


Finally  we  have 


|W2]x  - Ax| | 2 


. 86E-06 . 
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Conclusions 

As  far  as  accuracy  and  efficiency  are  concerned  the  generalized 
Lanczos-scheme  applied  to  skew-symmetric  matrices  is  more  attractive 
than  the  original  Lanczos-scheme  applied  to  the  squared  matrix,  with 
respect  to  both  accuracy  and  efficiency. 

T 

For  product-matrices  it  should  be  preferred  because  it  needs  no  LL  - 
decomposition.  However  special  care  has  to  be  taken  if  eigenvectors 
are  desired. 

One  of  the  main  difficulties  in  using  Lanczos-schemes  is  the 
monitoring  of  the  results.  This  difficulty  has  been  largely 
overcome  by  the  monitoring  process,  described  in  this  report. 

However  it  should  be  stressed  that  the  problem  of  determining  whether 
an  eigenvalue  of  the  original  matrix  is  single  or  not  has  not 
been  solved. 
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7.  Progrannature 


Subroutines,  in  Fortran  IV,  are  available  of  implementations  of  the 
Paige-style  Lanczos  and  of  the  generalized  Lanczos-schemes . These 
subroutines,  I.SVLAN  and  GENLAN,  as  well  as  the  subroutine  EVSCAN  are 
included  in  the  programl ibrary  ACCULIB  of  the  Academic  Computer 
Centre  Utrecht. 

In  this  section  documentation  and  listings  are  given.  This 


documentation  contains  a complete  example  of  use. 


TO  TRANSFORM  A SYMMETRIC  “ATPIX  A TO  TRIDIAGONAL  FORM  T,  BY  ORTHOGCNAU 
TRANSFORMATIONS  BY  THE  LA NCZOS-ME THO 3 . 

T ML  (ATP.IX  A H.t.OS  NOT  TO  BE  GIVEN  EXPLICITLY. 

EIGENVALUES  OF  T APPRCXI  'ATE  THE  EIGENVALUES  OF  A. 

AS  EACH  STEP  OF  THE  LANC73S  PROCESS  NEEDS  A MATRI X-VECTOR 
MULTIPLICATION,  THIS  FROCESS  IS  ONLY  SUITABLE  FOR  SPARSE  MATRICES. 


C HEADING  70  517 

c .................... 

C SUBROUTINE  LSVLANCN, IFIRST .M.fiESTRT , 01,  4X,SAVE0,U, CO, ALPHA. BETA) 

C LOGICAL  RESTRT 

C DIMENSION  Q1CN)  ,CJ (N) ,U(N>  , AL“HA (M) .BETA  CO 

C EXTERNAL  AX.SAVEC 

c 

C PURPOSE 

C 

c 
c 
c 
c 
c 
c 

c 

C INPUT-PARA lETERS 

C •••• 

C N -INTEGER.  THE  OFCER  OF  THE  FATRIX  A. 

C IFIRST  -INTEGER.  THE  FIRST  COLU tN  OF  THE  MATRIX  T,  WHICH  HAS  TC  BE 

C COMPUTED  CN  THIS  CALL  OF  LSVLAN. 

C ON  INITIAL  CALL  IFIRST  SHOULD  BE  1. 

C IF  LSV  LAN  IS  RESTARTED  (SEE  INPUT-PARAMETER  RESTRT), 

C IFIRST  SHOULD  BE  E 1 U A L TC  THE  LAST  COLUMN-NUMBER  OF  T IN 

C THE  PREVIOUS  CALL  = LUS  1. 

C H -INTEGER.  THt.  TOTAL  NUMBER  OF  COLUMNS  OF  T TO  BE  C0“PUTED 

C ( THE  NUMEL  R OF  COLUMNS  IN  EARLIER  CALLS  ARE  INCLUDED). 

C RESTRT  -LOGICAL. 

C RESTKTp. FALSE.  I INITIAL  CALL  FCR  A NEW  PROBLEM. 

C RESTRT-. TRUE.  1 INDICATES  A RESTART  AFTER  A PREVIOUS 

C CALL  CF  THE  SAME  PROBLEM. 

C Q1  -DIMENSION  Q 1 ( N ) • IF  RESTRT-.FALSE.  1 ARBITRARY  STARTING  VECTOR 

C FOR  THE  LANC70S  PROCESS,  NOT  NECESSARILY  OF  NORM  1. 

C IF  RESTRT =. TRUE.  « Q1(U)  SHOULD  HAVE  THE  SAME  VALUES  AS  CN 

C EXIT  OF  THE  PREVIOUS  CALL  (ALSO  CU T °UT -*>  A R AM ETER ) . 

C AX  -SUBROUTINE  AX(v,AY,N> 

C DIMENSION  Y (N) , AY (N) 

C THIS  USER-SUPPL IEB  SUBROUTINE  CELIVERS  FCR  A 

C GIVEN  VECTOR  Y THE  VECTOR  AY,  THAT  RESULTS  FROM 

C MULTIPLYING  THE  MATRIX  A WITH  THE  VECTOR  T. 

C Y SHOULD  NOT  BE  DESTROYED  WITHIN  Ax. 

C SAVEQ  -SUBROUTINE  SAVcQID.N) 

C DIMENSION  0(N> 

C THIS  USER-SUPPLIED  SUBROUTINE,  WHICH  CAN  BE  USED 

C TO  STORE  THE  COUU  INS  D CF  THE  CRTHOGONAU  TRANS- 

C FORMATION  MATRIX,  FOR  USE  IN  C C M DU  T I NG  THE 

C EIGENVECTORS  OF  A. 

C IF  EIGENVECTORS  ARE  NOT  CESIPEC,  THIS  SU"9CUTINE 

C MUST  STIUU  BE  S UpoL  IED-  -"l  T NEE  C NOT  ACTUAUUY 

C DO  ANYTHING. 


w 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 
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c 

c 
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c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


U -DIMENSION  JCN).  SCRATCH-ARRAY. 

QO  -OIMI-ENSION  DJ(N>.  IF  RES  TRT  = . r AL  SE  . « 01  NEED  NOT  TO  3E 

init iali2ed. 

IF  SESTET  :.  TRUE.  « OD  SHOULD  HAVE  THE  S A M £ VALUES  AS  CN 
EXIT  OF  THE  PREVIOUS  CALL  OF  LSVLAN 
( ALSO  CUT  PUT-PAR A 1 £ T E R ) . 

ALPHA  -DIMENSION  AL  °<&  (;i)  . IF  RE  S T RT  = . FALSE.  : NO  INITIALIZATION 

NECESSARY. 

IF  RESTR  T - , TRUE • « ALFHA(l)  UP  TC  AL°WA<MM>,  WHERE  WN 
1.1.1. LT.il)  IS  THE  VALUE  Oc  •!  IN  THE  PREVIOUS  CALL 
OF  LSVLAN  FOR  THE  S4m£  PROBLEM,  SHOULD  CCNTAIN 
THE  VALUES  OF  ALPHA  AT  EXIT  OF  THIS  PREVIOUS  CALL 
t ALSO  OUTPUT -PARAMETER) . 

ScTA  -DIMENSION  DETAIN).  IF  R£STRT=. FALSE.  J NO  INITIALIZATION 

NECESSARY . 

IF  RESTS!:.  TRUE.  « DETA(i)  UP  TO  SETA(;'M)  SHOULD  CONTAIN 
THE  VALUES  OP  BETA  AT  EXIT  OF  THE  PREVIOUS  CALL. 

( ALSO  OUTPUT-PARAMETER) . 


OUTPUT-PA-'.A  IETERS 


Q1  -DIMENSION  UIN).  U CONTAINS  THE  LAST  COLUMN  USEO  IN  THE  LANCZCS 

PROCESS  JF  THE  ORTHOGONAL  T RA  NSFO  R M AT  ION  MATRIX. 

( ALSO  INFUT-PARA-ET£R) . 

QO  -DIMENSION  DCIN).  Qj  CONTAINS  THE  COLU"N  PREVIOUS  TO  D1  IN  TH? 

transformation  process,  no  AND  Ql  are  necessary  for 
RESTART  PLR-’OSES  (ALSO  INPUT-PARAMETER)  • 

ALPHA  -DIMENSION  ALPHA  1. 1).  THE  DIAGONAL  OF  THE  TRIDIAC-ONAL  "ATPIX  T. 

WHICH  IS  ROUGHLY  SIMILAR  TO  A.  I ALSC  I NFUT  — PARAMETER)  . 

0 ET  A -DIMENSION  3ETAI,-).  THE  SU°  ERDIA  GONA  L ELEMENTS  CF  THE  MATRIX  T. 

3ETAC.I)  CONTAINS  RE  ST  ART  -INFORMATION. 

( ALSO  INFUT-PARA.-iETER)  . 


INTERNALLY  CALLED  SUBPROGRAMS 


VVIPP  (Z0  3Z3) 


REMARKS 


I)  LSVLAN  IS  SPECIALLY  OESIGNEC  FOR  THE  DETER  I IN  A T I ON  OF  THE 
EXTREME  EIGENVALUES  CF  A SPARSE  MATRIX  A.  IT  IS  OBSERVED, 

THAT  THE  EIGENVALUES  OF  TH  £ TRIO  I A GON  AL  M A T R I X T,  REPRESENTED 
BY  ALPHA  AND  3ETA.  PQ R INCREASING  ORDER  OF  T TEND  TO  THE  FIXED 
VALUES,  WHICH  CAN  DE  CONSIDERED  AS  EIGENVALUES  CF  A. 

FCR  PR03LEMS,  ARISING  WITH  7HE  DE  TE  R.'l  INA  T ION  OF  THE  EIGENVALUES 
WE  REFER  TO l 

VAN  KATS  J.M.,  VAN  OER  VORST  H.A., 

'•NUMERICAL  EXPERI  IENTS  OF  THE  PA  I Gr  -STYLE  LANCZOS  METHCO 
FOR  THE  COMPUTATION  OF  EXTREME  EIGENVALUES  OF  large 
SPARSE  MATRICES'*. 

1376,  ACCL,  TR  3 , 

II)  EIGENVALU_S  OF  THE  TRIDIAGONAL  MATRIX  T CAN  3E  COMPUTED  3Y 

IMTQLl 

ratdr 

BISECT 

ALL  CONTAINED  IN  EISFACK  t THIS  PACKAGE  IS  IN  UTRECHT  AVAILABLE 
AS  LIBRARY  ON  PERMANENT  FILES  L IN -'AS  ALGEBRA,  IO  = U ACCUI  ■ 

III)  LSVLAN  IS  SENT  IN  BY  J.  LEWIS  l COMPUTES  SCIENCE  DEPARTMENT) 
STANFORD,  U.S.A.. 


EXAKPLE  of  use 


IN  THE  FOLLOWING  PROCKAI  16  L ANCZOS- S T-"  PR  ARl  °ERcOP>ED  ON  A 71TH 
OROE S MATRIX  IWHICH  IS  KNOWN  AS  THE  MATRIX  WE J ♦ OF  WILKINSON). 
EIGENVALUES  OF  THE  MATRIX  W Cl  * ARE  COMPUTED  BY  EVSCAN,  USING 
THE  RESULTS  OF  LSVLAN. 


j 


1 
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PROGRAM  LNCZOS(CUTPUT) 

OI.-IENS  ION  TO  (21>,Q1(21),U(21>  » AL°HA  (16)  , BETA  (16  ) 
LOGICAL  FESTRT 
EXTERNAL  W21 =LS , 3AVEQ 
DIMENSION  TK(16,2) ,TK1(16,2> 

LOGICAL  LCW.UP.CIV 
N = 21 

IFIRST=i 
(1  = 16 

RESTRT  = . FALSE . 

NSEE 0= ..  1 0 7 77 
CALL  RAN ScT ( NS2ED) 

DO  1C  1 = 1, N 

10  Q1 (I) =RANF (SEED) 

CALL  LSVLAN(N,IFIAST,H,RESTRT,Ql,W21cLS,SAVEQ,U,Qi!. 

+ , ALPHA, BETA) 


"!T  = i 6 


CALL  E VSC  A N ( J , N , IT  , AL  PH  A , BE  T A , NEV  , LOW  ,l)P  , 0 IV  , 

♦ NCIV.TK.TKl, IERR.EFS) 

IF(IERR.NE.O)  STOP  ’’ERROR  IN  E IGE N VALUEC CNPU TAT  I ON ” 

I r { LOW ) PRINT  1000 
IP (UP)  PRINT  1310 
IF(DIV)  PRINT  1020, NDIV 
PRINT  1030, NEV 

C IF  NEV  EQUALS  3 NO  EIGENVALUE  INTERVAL  J HAVE  BEEN  OETECTED. 
C LSVLAN  SHDULO  EE  RESTARTED  FOR  FURTHER  DETECTION,  WITH  THE 
C DIMENSIONS  OF  ALPHA,  BETA,  TK  ANO  TK1  PROPERLV  ADJUSTED. 


IF  (NEV.E 

C.3)  GOTO  70 

DO  60  I=i 

,NEV 

60 

PRINT  13<,0,I,TK(I,H  ,TK(I.,2> 

70 

CONTINUE 

1000 

FORMAT (* 

COMVERGENCE  AT  LCWER 

SIDE 

») 

1010 

FORMAT (* 

CONVERGENCE  AT  UF°ER 

SIDE 

*.2) 

1020 

FORMAT (* 

*,I3,»  INTERVALS  AT 

LOWER 

END  «) 

1030 

FORMAT (* 

*,//,*  *,IJ,»  INTERVALS  ARE  FOUND  *,/, 

+ * 

N=>  L0WH330UN0 

UPPERBCUND*, /) 

IOmO 

FORMAT  (* 

13. ♦ *«2CE20.13«+ 

*)  ) 

END 


SUBROUTINE  SAV£C(Q,N) 

DIMENSION  Q ( N) 

C NO  EIGENVECTORS  A RE  COMPUTED,  SO  WE  DO  NOT  STORE 
C THE  ORTHOGONAL  TRANSFOR-lAT  ION-MATRIX. 

RETURN 

END 


SUBROUTINE  W 21 PLS ( X , A X , N ) 
DIMENSION  X(N) ,AX(N> 


C MATRIX  IS  TAKEN  FROM  WILKINSON  (THE  ALGEORAIC  E IGENVALU  E PRO  EL  EM  ) 
C ANO  IS  CALLED  W21*. 


C 

C 

C 

C 

C 

C 

c 

C 


10 

0 


i o 

9 10 

i a i a 


18  10 

0 19  1 

0 110 


,C  best  QU^iri  lfRACI1“Aj 
THIS  HAOE  IS  — 10  DDfi  — 1 ^ 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
A X (H  =10  • X ( 1 • + X (2 ) raOM  COPY  FURNISHED  TO  DLC  , 

00  10  1=2,23 

AX(I>  = X(I-l)+I4PS(ii-I>*Xm*X(I*l» 

10  CONTINUZ 

4 X (2  i ) =X  (2j)  ♦10*XC21) 

RETURN 

END 

THE  OUTPUT  OF  THIS  PROGRAM  IS* 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CONVERGENCE  AT  LCWER  SIDE 
CONVERGENCE  AT  UPPER  SIDS 

2 INTERVALS  AT  LOWER  END 


6 INTERVALS  ARE  FOUNC  . 
NR  L0WER30UN0 

1 -.11 25rL1E221-9E»3i 

2 .25330581709739+03 

3 . 7303952209h22E+G1 

A . 8 33 39*11 223762+01 

5 . 9213678e<.7337E*ol 

6 . l07RbigRjL8290E*32 


UPPER30UND 

-. 112544l5221l9E*0l 
.2533J53173996S+00 
.'7C0395220942‘.E+G1 
. 3 j3394ii223765+ 01 
.9213578647337E+01 
.10746i94lS290S+02 


MET  HO  0 


BASED  ON  THE  PAIGE  STYLE  LANCZOS  PROCESS  DESCRIBED  IM 
LEWIS  J.,  THESIS  (TO  APPEAR!. 

VAN  KA73  J.N.,  VAN  DER  VORST  H.A.,  "NUMERICAL  EXPERIMENTS  OF 
THE  PAIGt-STYLE  LANCZOS  FETHCC  F C R THE  COMPUTATION  CF 
EXTREME  EIGENVALUES  CF  LARGE  SPARSE  MATRICES", 

1976,  ACCU  TR3. 

VAN  KATS  J.  I.,  VAN  OER  VORST  H.A.,  "AUTOMATIC  MONITORING  CF 
LANCZOE  cCME.lES  FOR  SY'-'METRIC  AMO  SKEW-SYMMETRIC 
GENERALIZED  EIGENVALUE  PRC3LEHS,  19ZZ,  ACCU  TR7 . 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  LSVLAMN,  FIRST,  M,  RESTRT,  Qi,  AX,  SAVED,  1J,  DO , 

1 ALPHA,  BETA! 

SUBROUTINE  TC  CARRY  OIJT  LANCZOS  PROCESS  FOR  A SYI:, METRIC  REAL  MATRIX 
INR'JTt  N =>  THE  SIZE  OF  THE  MATRIX 

FIRST  =>  THE  FIRST  COLUMN  TC  BE  CO’-RUTEC  0 N THIS  CALL  to  THE 
THE  SUBROUTINE  (ON  INITIAL  CALL,  FIRST  WILL  BE  1) 

I =>  THE  I A X IilU. N Nlf-EZR  OF  COLLINS  CF  THE  PROCESS  TC  BE 
COMPUTED  • 

EES  T R T =>  A LOGICAL  VARIABLE  incicating  whether  we  are 
RESTARTING  FRCI  A PREVIOUS  PARTIAL 
T RIDIAGOMLIZAT  ION  or  STARTING  a NEW  ONE. 

Di  =>  THE  FIsST  COLUMN  OP  1,  A REAL  VECTOR,  NOT  ASSUME 
TO  BE  CF  NCRf  1. 

AX  •=>  A S J TROUT  I Nc  TO  CARRY  OUT  THE  J A T R I X - VEC  TOR 

MULTIPLICATION.  INPUT  TO  AX  IS  A REAL  VECTOR. 

OUTPUT  SHOULC  BE  A PEAL  VECTOR. 

S A >/  E Q =>  A SUBROUTINE  WHICH  SAVES  THE  COLUINS  GENERATEC  BY 
THE  LANCZCS  ALGORITHM,  FOR  USE  IN  COMPUTING 
EIGENVECTORS.  IF  EIGENVECTORS  ARE  NOT  DESIRED, 

THIS  SUBROUTINE  MIST  STILL  BE  SUPcLIEn--IT  NEED 
19  T ACTUALLY  DC  ANYTHING. 
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SCRATCH:  U =>  A real  \JZ CTOk  Cf  LENGTH  N 

G.j  =>  A REAL  VrCTCR  0^  LENGTH  N 

NEITHER  IJ  M3  03  N££  G n~  I N I T IAL  I 7£ 0 
TH^  TO  ITlNTS  OP  3l,  )0  AMI  U WILL  EE  DESTROYED 

OUTPUT*  AL  °HA  -•»  THE  HAGCHAL  C ~ THE  TR 10 1 AGONAL  MATRIX  Wt-ICH  IS 

ROUGHLY  SI  I L A < TO  A. 

3ZTA  =>  THE  S’.MH  ACIAGGNAL  '3c  THE  TR  ID  I AG  C N AL  MATRIX 
ALPHA  AM  GET  A ARE  PEAL  VECTORS. 


PRUC.-.S3  IS.C3*  T h I 3 SU  GROUT  I M I “ ° L £ T I NT  S A R-~FINEO  VERSION  CP  THE 
-’A  IG  1 - G T x LE  LA  NC  US  PROCESS. 


I \ T _ r,  I R \ , ,•  f FIRST 

L Ot> C TL  ~d  3 T U 
EXTERNAL  ax,  SAVI3 
real  U(u>.  TiCI) 

3 w A L J < N > 

R : A L ALOHA  ( i)  , GET  A ( i) 


REAL  LTiE.)F»  LALPHA,  LOETA 
REAL  E ° 1L  CN 
1 1\  T EGER  I,  < 


IACHINE  DEPENDENT  3UANT IT  IES  * 

Z^SLON  13  THE  RELATIVE  T RU  NC  A T ICN  LEVEL  CF  SINGLE  PRECISION 
IMA  . ! n 3 L C N / A , E - L 4 / 


IF  (RESTRT)  GC  TO  60 


I IITIALI7ATION 


STEF  i?  NORlALIZE  1i  TC  UNIT  LENGTH 
L Tc . ; P = VVIPP(  U,i,:U,l,N,T) 

IP  f LTEOP.EQ.^.  ,1)  GO  TC  £U 
L TE . iP  = i. J/SQRT  ( LTEMP) 

00  _ J I =1,0 

u ( I)  = Ql  d>  *LT£.-1P 
1J  CONTI  NUC! 

C CO  I CUTE  Q2 


C 

c 

c 


c 


EJ  CALL  GAVECn.,  N) 

CALL  A X ( U , U,  N) 

L AL  ->HA  = VVIpP (U , i , Ql , i * N , T ) 
Cl  0 3 J I = i , N 


U(I)  = IJ(  I)  - L A L PH  A * C L ( I» 

TO  CONTINUE 

FOR  USE  III  AN  ITERATIVE  LANCZOS  PRCCESS,  SINCE  GETA  ^flY  GZ 
VERY  3 lALL,  WE  RE C RTHOG ON AL I ZE  02  WITH  RESPECT  TC  01  IN  ORDER 
TO  AVHJ  LOSS  OF  J LCT  OF  CLS  ORTHOGONAL  ITY  RIGHT  OFF. 

L TE,  IP  = VVIPPUJ,  i,OI,  1,  N ,T> 

0 0 t 0 1=1,3 

HI)  = U C I ) - LT.Efp’Cid) 

'♦j  CONTINUE 
NOW  NORlALIZE  02 


L 3E  T A = VVIPP(U,i,U,i,N  , f) 
LOETA  = SORT  T LOETA) 

L TE  IP  = i.O/LOET  A 
NO  E>  J 1 = 1,1 

OJ(I)  = OKI) 

Oi(I)  = U (I) * L T £ ! I P 
5 J C ON  T I NtJiE 

CALL  SAVES  ( II , N) 


THIS  TAGE  IS  to  UlKl  _ — ^ 


1 
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0 


c j i »ijTc.  i<itial  a ’in  oeia. 

A w ° i < A (I)  = L A L.  A 

11TA  (_)  - L.  1-TA 

t<  - 2 
3 J r 0 7 j 


THIS  i'ACK  IS  HS.  T Q"M  I ?Y  RA  ' 
Oil  \ KUKMbitu.  lo  OJJ 


ncAivu 


c ============================================ 

CRL3TARTING 

c ============================================ 

C THE.  u.=  2 fJST  HAVE  p-.STCFJO  31  AMC  C3,  AL^UA.  1'TA  ANO 

C nA»^A  A3  THEY  :(£•*_  AT  THE  END  CF  TH?  LAST  CALL.  L3=TA  IS 

C ST  ) O IN  3ETA (PIRST-.I 

o»  < - Ft?ST 

L.3ETA  = 3 E T A ( K - * ) 

CALL  SAVEQ(U»  N> 

C 

C NOW  CA  NY  U)T  LANC20S  "’ROC-JSS  FOR  < = 1,  h,  ....  HJ 

c 

C F )R  < = 2 TO  --1 

C < CO  : 11 L T E SlICCc.LUING  COHJ  - IS  OF  ORTHOGONAL  “ATSIX  0 > 

• o iF  ( X , 1 » T , i%  ) (.» C (0  1 0 ij 

C COIFUTE  °'fVICU3  AL°HA  AIIJ  CIEECTICJ  CF  N^XT  VECTOR 
LTK  IP  = VVPPI  U,i,  U,i,N,T) 

CALL  AX(  lx,  U,  fT ) 

LAL-'LHA  = VVIP=(II,1,11,1,N,T» 

A LPH  A ( K ) = L A L aH  A 

TO  T J 1 = 1,  N 

JO  = U ( I » - LALPHA*GL  < I > - L3ETA*Q3(U 

.33  CJNTI  IDE 

C NOXrALIZ?  NEXT  V:.CT0R 

LOETA  = VVIPF  (II,  1,11,1  ,N  ,T» 

Lt)c  TA  = SQ  ET  (L3ITA) 

IF  (L3£TA.LT._-’SLCN)  LPSTA  = ZFSLGN 
LT£.|P  = i.U/LOETA 
OlT  A (<)  = LHE  T A 

TO  lj  I = * , N 

U ( I)  = OKI) 

11  ( I)  = (J  ( I)  *LT£  IP 
HJ  continue 

IF  (K.LT.M  CALL  SAVtQCCl,  N) 

K = K * 1 
GO  TC  70 
1 ) j RETURN 
E NO 


C MEAOING  7U51S 
c 

C SUBROUTINE  GENLANCN,  I F I RS T , , RE  ST  RT  , 0 1 , CX  , OX  , S AV E 3 , U ,U  W , 

C QD,  ALPHA. EETA.AKTU 

C LOGICAL  RESTRT.APTI 

C DIMENSION  Ul  C N»  ,Cj  (N)  ,U(NI  , ALPHA  (.!>,  BETA  (B)  ,UW(N) 

C EXTERNAL  CX.DX, SAVED 

c 

C PURPOSE 
c 

C TO  GENERATE  A TRI DIAGONAL  MATRIX  T OY  THE  LANC70S  PROCESS  FOR  A 

C MATRIX  C»0,  WHERE  C IS  WHETHER  SYMMETRIC  OR  SKEW- SYMMETRIC  A KD 

C 0 IS  SYMMETRIC  PCSIT1V--  DEFINITE, 

C THE  MATRICES  ‘1  AND  C MED  NOT  TO  Oc'  GIVEN  EXPLICITLY. 

c eigenvalues  of  t appkcxi  iate  the  eigenvalues  OF  C*B. 

C AS  EACH  ST  >’  OF  th;  lAHCZOS  PRCCts;  NEEDS  2 H a TR  i x-Vt  C TOR 

C MULTIPLICATIONS,  TMIS  PRJCESS  IS  V_RY  ATTRACTIVE  FOR  SPARSE 

C MATRICES. 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


INPUT-PARAMETERS 


N - INT  EGER . THE  ORCER  OP  THE  MATRICES  C AND  B. 

IFIfiST  -INTEGER.  THE  INCEX  OF  THE  FIRST  COLUMN  OF  THE  MATRIX  T,  WHICH 

HAS  TO  HE  CO..PJTE0  ON  THIS  CALL  OF  GENLAN. 

ON  INITIAL  CALL  IFIRST  SHO'JLO  3c  1. 

IF  GENLAN  IS  RESTARTED  (SEE  I N PU T - P A R A ME T E S RESTRT), 

IFIRST  SHOULD  3E  EQUAL  TO  THE  LAST  COLUMN-NUMBER  OF  T IN 
THE  PREVIOUS  CALL  PLUS  i. 

M -INTEGER.  THE  TOTAL  NUMBER  OF  C 0 LIMNS  OF  T TO  BE  COMPUTEO 

( THE  NUMEER  OF  COLUFNS  IN  EARLIER  CALLS  ARE  INCLUCE'Dt  . 
RESTRT  -LOGICAL. 

RESTRT =. FALSE.  « INITIAL  CALL  FOR  A NEW  PROBLEM . 

RESTRT=. TRUE.  » INDICATES  A RESTART  AFTER  A PREVIOUS 
CALL  CF  THE  SANE  PRO 3LE f* . 

Ql  -DIMENSION  Ql(N>.  IF  R£STRT=. FALSE.  « ARBITRARY  NCN7ERO  STARTING 

VECTOR  FCR  THE  LANCZOS  PROCESS,  NOT  NECESSARILY  CF  NC Rx  1. 
IF  RESTRT=. TRUE.  x Qi  (N)  SHOULD  HAVE  THE  SAHE  VALUES  AS  CN 
EXIT  OF  THE  PREVIOUS  CALL  (ALSO  OUTPUT-PARAMETER!  . 

CX  -SUBROUTINE  CX(Y,CY,N) 

DIHENS  ION  Y ( N ) ,CY(N) 

THIS  USER-SUPPLIED  SUBROUTINE  CELIVERS  FCR  A 
GIVEN  VECTOR  Y THE  VECTC.R  CY,  THAT  RESULTS  FROx 
THE  MULTIPLICATION  OF  Y BY  THE  (SYMMETRIC  OR 
SKEW-SYMMETRIC)  MATRIX  C. 

Y SHOULD  NCT  BE  OESTROYEO  WITHIN  CX. 

BX  -SUOROUTINE  BX(Y,3Y,N) 

DIMENSION  Y(N),3Y(N) 

THIS  US£R-SU00L IEO  SUBROUTINE  CELIVERS  fcp  A 
GIVEN  VECTOR  Y THE  VECTOR  BY,  THAT  RESULTS  FROx 
THE  MULTIPLICATION  OF  Y BY.  THE  SYMMETRIC 


™IS  PAG?  IS  BEST  QUALITY  PRACTICABLE 

KAW  COPY  FURNlSHjl:  TO  l DC 


C 

c 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C 

c 

c 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


POSITIVE  DEFINITE  MATRIX  B. 

Y SHOULD  NOT  BE  DESTROYED  WITHIN  3X. 

SAVEQ  -SUOROUTINE  SAVEQ(Q.N) 

0 1.IE  NS  ION  Q ( N) 

THIS  USER-SUPPLIED  SUBROUTINE,  WHICH  CAN  BE  L'SEO 
TO  STORE  THE  COLU1NS  D OF  THE  CRTHCGONAL  TRANS- 
FORMATION I.ATRIX,  FOR  USE  IN  COMPUTING  THE 
EIGENVECTORS  OF  C (IN  THIS  CASE  0 SHOULO 
BE  THE  IDENTITY  MATRIX). 

IF  EIGENVECTORS  ARE  NOT  OESIRED  CR  WHEN  D IS 
NOT  THE  10 ENTITY  MATRIX,  THIS  SUBROUTINE 
MUST  STILL  BE  SUPPLIED — IT  N EEC  NOT  ACTUALLY 
DO  ANYTHING. 

U -DIMENSION  U(N> . SCRATCH-ARRAY. 

UW  -DIMENSION  UW(N>.  SCR A TCH- AR R AY. 

QQ  -DIMENSION  QV(N).  IF  -RSSTRT  = . FALSE  . « Ql)  NEED  NOT  TO  BE 

INITIALIZED. 

IF  RESTRT:. TRUE.  « DO  SHOULD  HAVE  THE  SAXE  VALUES  AS  CN 
EXIT  OF  THE  PREVIOUS  CALL  OF  GENLAN 
( ALSO  OUTPUT-PARA  IETER) . 

ALPHA  -DIMENSION  ALPHA (M) . IF  RE S T R T = . F A L SE . » NO  INITIALIZATION 

NECESSARY. 

IF  RESTRT=. TRUE.  X ALPHA ( 1 ) UP  TO  ALPHA(XH),  WHERE  MM 
(NM.LT.M)  IS  THE  VALUE  OF  M IN  THE  PREVIOUS  CALL 
OF  GENLAN  FOR  THE  SAXE  PROBLEM,  SHOULD  CONTAIN 
THE  VALUES  OF  ALPHA  AT  EXIT  OF  THIS  PREVIOUS  CALL 
( ALSO  OUTPUT-PARA  DETER) . 

BETA  -DIMENSION  BETA (M . IF  RES TRT=. FALSE . 1 NO  INITIALIZATION 

NECESSARY. 

IF  RESTRT:. TRUE,  j 8 £ T A ( * ) UP  TO  GETA(MM)  SHOULD  CONTAIN 
THE  VALUES  OF  BETA  AT  EXIT  OF  THE  PREVIOUS  CALL. 

( ALSO  OUTPUT-PARAMETER) . 

ANTI  -LOGICAL. IF  THE  MATRIX  C IS  SKE W -S YM MET R IC,  ANTI  SHOULD  BE  SET 

TO  .TRUE.,  IF  THE  MATRIX  C IS  SYMMETRIC  ANTI  SHOULO  BE 
SET  TO  .FALSE.. 


OUTPUT-PARAMETERS 


Ql  -DIMENSION  D1(N).  Ql  CONTAINS  THE  LAST  COLUMN  USEO  IN  THE  LANCZCS 

PROCESS  CF  THE  ORTHOGONAL  TRANSFORMATION  MATRIX. 

( ALSO  INFUT-PARAMcTER) . 


THIS  PAGE  IS  BEST  QUALITY  I’RACI'iCA&Ull 
KKvM  Gv)i  Y EUKhlSiUL*  iy  L\DO  - 


c QO  -OIHENSION  Tl  (N)  . f)J  CONTAINS  THE  COLUMN  PREVIOUS  TO  01  IN  THE 

c TRANSFORMATION  FRGOESS.  TO  AND  01  ARE  NECESSARY  FOP 

C RESTART  PURPOSE  S (ALSO  I N 3ll  T- PAR  AF  £ T£  ? t . 

C ALPHA  -DIMENSION  ALPHA!  O.  THE  DIAGONAL  OF  THE  TSI DIAGONAL  MATRIX  T, 

C WHICH  IS  ROUGHLY  SIMILAR  TO  C*B.  ( ALSO  IN  °U  T-P  A c A .*•  £ T E S ) . 

C BETA  -DIMENSION  3ETAO).  THE  Sll  P E °D  I A G DNA  L ELEENTS  CF  THE  HATCI*  T. 

c 9E  T A ( M ) CONTAINS  RE  START-IN FORMATION. 

c IF  C IS  SKEW-SYr  1ETAIC,  THEN  THE  SUDOIAGCNAL  ELEFENTS 

C ARE  EQUAL  OUT  OPPOSITE  IN  SIGN. 

C ( ALSO  INFUT-PARA.iETERI  . 

c 

C INTERNALLY  CALLED  SUBPROGRAMS 

c .................... 

C VVIPP  (70929) 

c .................... 

C REMARKS 

c .................... 

Cl)  IF  0 IS  THE  IDENTITYtlATRIX  AND  C IS  SYMMETRIC,  THEN  LSVLAN 

C (7U517)  SHOULD  EE  PREFERRED  BOTH  FOR  TI.-E  AND  STORAGE  CONSIOEKA- 

C TIONS. 

C II)  GENLAN  IS  SPECIALLY  ATTRACTIVE  FOR  THE  DETERMINATION  OF  THE 

C EXTREME  EIGENVALUES  OP  A SPARSE  FA  T R IX  C*B  WITH  C AND/Oc  n SPARSE. 

C THE  EIGENVALUES  CF  THE  TRIOIAGONAL  MATRIX  T.  REPRESEN TEC 

C BY  ALPHA  AND  BETA,  FOR  INCREASING  ORDER  OF  T TEND  TO  THE  FIXEC 

C VALUES,  WHICH  can  BE  CONSIDERED  AS  EIGENVALUES  CF  C*B. 

C FOR  PROBLEMS,  ARISING  WITH  TM:.  DETERMINATION  OF  THE  EIGENVALUES 

C HE  REFER  TO) 

C VAN  KATS  J.H.,  VAN  0”  R VORST  M.A., 

c “numerical  experiments  of  tm_  paig; -style  lanczos  f e t h c b 

C FOR  THE  COMPUTATION  OF  EXTREME  EIGENVALUES  OF  LARGE 

C SPAP.SE  FAT  RICES", 

C 1976,  ACCU,  TR5. 

C III)  THOSE  EIGENVALUES  OF  THE  TRIOIAGONAL  MATRIX  T,  REPRESENTED 

C BY  THE  ARRAYS  ALPHA  AND  BETA,  WHICH  CORRESPOND  TO  EIGENVALUES 

C OF  THE  ORIGINAL  FROOUC T-.IA TR IX  C*-3  CAN  BE  COMPUTED  BY  EVSCAN 

C (70519). 

C IV)  IF  THE  MATRIX  IS  SKEW-SYMMETRIC,  THEN  THE  EIGENVALUES  DETERMINED 

C BY  EVSCAN  (7U51E)  SHOULD  BE  MULTIPLIED  BY  SORT(-l). 

C V)  GENLAN  IS  WRITTEN  BY  H.A.  VAN  9ER  VORST  (ACCU)  UTRECHT. 

C .................... 

C EXAMPLE  OF  USE 

c .................... 

C 

C THE  FOLLOWING  PROGRA:  COMPUTES,  WITH  EVSCAN  (70519),  THE  EIGENVALUES 

C OF  A SKEW  SYMMETRIC  CAT.RIX. 

C 

c 

C PROGRAM  LNCGEM(OUTPUT) 

C DIMENSION  Cu  ( 20  ) ,91  ( 2 J ) ,U  ( 20  > ,UN  ( 2 0 ) * 

C * ALPHA  (25)  .BETA (25) 

C LOGICAL  R2STRT, ANTI, LOW, U°.DIV 

C EXTERNAL  1 K E W A , ICENT I, SAVED 

C DIMENSION  TK (25,2) ,TK1 (E5,2) 

C N = 2U 

C IF  I RS T = 1 

C M=25 

C RESTRT=. FALSE. 

C ANT  I = • T RUE , 

C NSEEO-JC 0 E77 

C CALL  RANSET (N3EEC) 

C DO  ID  I=i, N 

C 10  q„  (I » =RAN f (SEED) 

C C A NON7ERO  RANCC.I  STARTVECTQP 

C CALL  GtMLAN(N,IcIRST,F,  RE  STRT,  (91, SKEWA.IOENTI, 

C • SAVED, U,UW,Q0, ALPHA, D.TA, ANTI) 

C MT =25 

C CALL  EVSCAN(  l,N,FT,ALFHA,BETA,NEV,LOW,LIP,riV, 

C • NO  IV, TK , TK1 , I£RR«E°S) 


r 
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C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IF  ( I ERR.  He.  jl  STOP  "ERROR  IN  EIGENVALUE  COMPUTATION" 

IF  (LON)  FAINT  1300 
IF  (UP)  PRINT  1010 
IF  (OIV)  FAINT  1020* Nn IV 
PRINT  10  30  , NEV 

C IF  NEV  EQUALS  0.  NO  EIGENVALUE  INTERVALS  HAVE  SEEN  CETeCTEG, 
C GENLAN  SHOULC  RE  RESTARTED  FOR  FURTHER  INFORMATION,  WITH  THE 
C 0 1 ME  NS  I ON S OF  ALPHA, OETA,T<  ANO  TK1  PRO°ERLV  ADJUSTED. 

IF  (NEV.EQ.u)  GOTO  70 
OO  60  I -i « NEV 

60  PRINT  lut,  ,I,TK(I,i) ,TK(I,2) 

70  CONTINUE 

1000  FORMAT!*  CONVERGENCE  AT  LOWER  END  ») 

1013  F CRilAT  { * CONVERGENCE  AT  U = PER  ENC  *) 

1320  FORMAT ( * * , I 3 , * INTEPVALS  AT  LOWER  END  *) 

1030  FOR.lAT(»  *,//,*  »,IJ,*  INTERVALS  ARE  FOUNC  *,/, 

• * NR  LOWERROUND  U F PE  RBOU  NO  *,/) 

1040  FORMAT ( * * , I 3 , * *,2(F20.13,*  *)) 

END 


SUBROUTINE  SAVE 3 (0 » N) 

DIMENSION  C ('!) 

C NO  EIGEfIVECTCRS  ARE  REQUIRED,  SO  WE  DO  NOT  STCRE  THE 
C ORTHOGONAL  TR A NSFORH A T I ON  MATRIX. 

RETURN 

ENO 


SUBROUTINE  SKEWA  (X,CX,N) 

OIUENS ION  X C M ) ,CX(N) 

C THE  SUBOIAGCNAL  ELEMENTS  OF  THIS  SKEW  SYMMETRIC  l,ATRIX  ARE 
C CHOSEN  TO  9£  -1,  THE  SUPE  RO I AG  ON  AL  ELEMENTS  ARE  «■!  ANC  ALL 
C OTHER  ELEMEMS  ARE  0. 

CX  (1)  = X (2) 

Nl=N-i 

DO  10  1=2, M 

10  CX  (I)  =X  ( I <1) -X  (I-l> 

CX (N) =- X ( Ni  ) 

RETURN 

END 


SUBROUTINE  IDENT I ( X, 9X , M) 

DIMENSION  X(N),;iX(N) 

C THE  TRANSFORMATION  9 IS  THE  IDENTITY  “ATRIX  IN  THIS  CASE. 

OO  ,U  1=1, N 
10  9X(I)=X(I) 

RETURN 

ENO 


THE  OUTPUT  OF  THIS  PROGRAM  ISI 


CONVERGENCE  AT  LCWER  ENO 
CONVERGENCE  AT  UFPER  ENG 
20  INTERVALS  AT  LOWER  ENO 


t 

CO!  T FUKMsHJio  IO  DDC 


20 

NR 

INTERVALS  ARE  FOUND 
LOWEROOUNC 

UFPERBCUND 

1 

-1.9776616524499 

-1.9776616524498 

2 

-1. 3111456115717 

-1.9111456115716 

3 

-1. 8019377355145 

-1.8019377358044 

4 

-1.6524775436317 

-1 . 6524775436316 

5 

-1.4661J374 36591 

-1.4661037436593 

6 

-i .2469756037172 

-1.24697  96  0371*2? 

7 

-. 9999993339350 

-.9999939599999 

8 

-.  73  J632  J-»  37327 

- .730682J-.  37327 

9 

-. 4450413679125 

-.4450419679125 

10 

-. 1494601371729 

-.1 494631 971728 

noooooooooooooooooo 


11 

.1491*60187172.1 

.1494601 971723 

12 

, 4*,  5 J 4 1 3o7  91  2 5 

.4453416679125 

13 

,7306320417326 

.7336820487326 

m 

. 999995H9  33 »9o 

.9939999999  396 

15 

1.24697963  37173 

1.2469796037173 

16 

1. 46610  27436592 

1.4661037436593 

17 

1. 6524775436313 

1.65247754  36  314 

18 

1,88^.9377358340 

1. 8019377358  343 

19 

1.9111456115715 

1.9111456115716 

20 

1. 9776616524493 

1.9776616524493 

NET HOO 


THE  GE  N=RA  LIZSO  LANCZC3  SCHEIE  AS  DESCRIBED  IN* 

VAN  KATS  J.N.,  VAN  OER  VORST  M.A., 

"AUTOMATICAL  ,-!OM  TOR  ING  OF  GENERALIZED  SYMMETRIC 
OR  SKEWSV ('METRIC  IANCICS  SCHEMES", 

1977,  A CCU ■ UTRECHT,  TR7. 


ALPH A (K»  , OETA (M  > 


AX 
0 R 


3UR  ROUTINE  GE.mLANO,  IPIAST,  I 
1 .'j  » ALPHA,  rl_:TA,  ANTI) 

LOGICAL  AMI,  RESTRT 
EXTERNAL  AX,  EX,  3AVE0 
REAL  Ql(N),  Oj('I),  U(N>,  UW(N) 

REAL  LTf  IP,  L)Ln,U,  L OS  T 0 
REAL  P3LC  I 
I N T . G i R I,  K 
DATA  EPSL  CN  /i.E-14/ 

GEN  LAN  GENERATES  RO.I3  CCR.  THE  PROOICT  MATRIX  3X 
OX  IS  SY.1IETRIC  A NO  POSITIVE  DEFINITE, 

IS  WHETHER  SY.-h£TR  1C  t AMT  I = . FALSE  . , 

ANTI'SY  lUETSIC  J ANTI^.TRLE. 

IF  (R:_STRT)  GO  TO  70 
CALL  OX (01,  UN,  N) 

LTE:IF  = VVIPP  (-’Uf  i»'.)W  ,1  ,N,T) 

IF  ( LTE  IP.EQ, .. J GO  TO  EJ 
LTE.'P  = 1.U/S0RT  (LTE.HP) 

00  _J  1=1, N 

U(I)  = Q1(I)  *LTI,|P 
UW(I)  = UW(  I)  *LTE  T° 

CONTINUE 

i CALL  SAVEQCQl,  N) 

CALL  mX(UW,  li,  N> 

IF  (AUTI)  GO  TO  hv 
L AL:,HA  = VVIPF(U,i,UW,l,N,T) 

TO  JJ  1 = 1 ,rj 

U ( I)  = U(I)  - LALPHA*Gl(I) 

CONTINUE 
CONTINUE 

CALL  3X(U,  IJU,  Nl 
L T E . i D = VVI°P  00,1,01, 1,N,T) 

00  '» u 1 = 1 )N 

UCD  = (1(1)  - ITE,|P*C1  ( I) 

CONTINUE 

LRETA  = VV  IPP  ON  , 1 ,11  , 1 , N ,T) 

LOETA  = SORT (10  ETA) 

L T E i 3 = 1.0/LOE7A 
00  dJ  1=1, N 

0 J ( I ' = 01(1) 

Old)  = IJ  ( I ) -LTE, IP 
JH  ( I)  = UW(I)*LT.E'1P 
b'J  CONTINUE 


RE  ST  RT , Cl,  AX,  S X , SAVED,  U,  UL, 


AX,  LHER? 


JJ 

j 


6J 


<0 


cc 


V 

N>v 


9 vO' 


C-LL  '.AVEQCU,  )» 

IF  ( A J T I » LAL'UA  = i. 

ALPHACU  = LAL  ’HA 
T A ( _ ) = L3ITA 

< = J 
G 0 TO 

rj  < = ifir-it 

LIETA  = Bt  T A (K-_  ) 

caul  ;av£CUi,  n 

j J IF  (K.GT.i*)  GC  TO  1 m J 

call  ..xcjw,  u,  n 

IF  CA'ITI)  GO  TO  1J.J 
L AL  3MA  = WIFF(lJW,l,UtitNfT> 

ALFO  A ( < ) = L AL  PH  A 

TO  lj  1=1 , 0 

t)(I>  = U(T)  - LALPH  4*Ci  ( I ) - LBFT  A*Q  3(1) 


)J 

CONT 

IN  UL 

go  r 

0 ICO 

13  j 

Al^HACK)  = 

U • 

no  _ 

A.  U I - I 

, N 

• ) c r ) = 

'J  ( I ) f LGE 

1-3 

CONT 

I NUF 

12J 

call 

3XCU, 

JN,  N ) 

L IF  T A = VVIPPUJ,  1,UH,1,  ► .T) 

LOETA  = SQUCLTETA) 

IF  (LOIT A.LT. IPGLOrO  L 3 c 7 A = I F 3 L Q N 


S Jfp®  IS  REST  QUAim  ^iACn  CABLE 
■WJOiM  COFI  io  QDC 


30 


1 m0 


LTF.  |P  = 1 . J/L3£T A 
RET  UK)  = L3FTA 
3 0 .j  i I = _ , N 

l j ( i ) = q i m 
U<  I>  = U < I)  *LT£.I° 

)W(I)  = UNCI)  *L  T£:‘l  3 

CONTINUE 

IF  CK.LT.>)  CALL  SAVEQCGl,  N) 
K = X *■  1 
GO  TO  3J. 

RGTiJR  ) 

1 NO 


c .................. 

C HEADING  7051 9 
c .................. 

C SUOROUTINE  £ VSCA  N ( .'1 , N ,MT  . AL  PH  A , 3cT  A , N£V,  LOW,  UP  , 0 IV,  NO  I V , TK, 

C TK1 , IERR.EPS) 

C DIMENSION  ALPHA  CM  , DET  A (H)  ,TK(HT  ,2)  »TK1(MT,2I 

C LOGICAL  LOW,UP,0IV 

c .................. 

C PURPOSE 

c .................. 

C TO  DES  T ILL  AT  I A ROW  CF  DISJUNCT  INTERVALS,  EACH  OF  WHICH  CONTAINS  AT 

C LEAST  ONE  EIGENVALUE  CF  A GIVEN  TR 10  I AGON AL MA T RI X T,  WHICH  ARISES 

C IN  THE  SINGLE-VECTOR  LANCZOS  PROCESS. 

C 

C INPUT-PARAMETERS 

c .................. 

C M -INTEGER.  THE  ORCER  OF  THE  TSIDI AGONAL  MATRIX  T.  AS  GENEP.ATEC 

C BY  THE  LANCZOS- PROCESS. 

C CM  IS  EQUAL  TO  THE  NUMBER  OF  ITERATIONS  IN  THIS  FRCCESSI, 

C N -INTEGER.  THE  ORCER  OF  THE  CSTGINAL  MATRIX  A,  ON  WHICH  THE 

C LANCZOS-FSOCESS  MAS  EEEN  APPLIED. 

C MT  -INTEGER.  THE  NUMBER  OF  ROWS  IN  THE  ACTUAL  DECLARATION  CF  THE 

c ARRAYS  TK  AND  T Ki • 

C ALPHA  -DIMENSION  ALPHA  «,!)  . THE  DIAGONAL  OF  THE  TRIOIAGONAL  MATRIX  T 

C AS  OELIVEREO  3Y  THE  LANCZCS-PROCESS.  THE  ELEMENTS 

c OF  ALPHA  ARE  NOT  ALTERED. 
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-DIME  NS  ION  .Di_TA(fl.  QETA(i)  THROUGH  TJTaC--l.»  CONTAIN  TH£ 
SUPERJ  IAGCNAL  ELEMENTS  OF  THE  TRIDIAGONAL  “ATRXX  T. 
THE  ELEMENTS  OF  RET  A ARE  NOT  ALTERED. 

-OIMENSI  ON  TK.l(.*>r,2).  A SCRATCH-ARRAY. 


C OUT  FU  T- P A i A lETcRS 

c .................. 


-INTEGER.  THL  NUJ-RER  OF  DIFFERENT  EIGENVALUE  INTERVALS,  IK 

WHICH  EIGENVALUES  OF  THE  ORIGINAL  MATRIX  A ARE  CONTAINED. 
-LOGICAL.  IF  LOW  IS  .TRUE.  » CONVERGENCE  AT  THE  LOWER  SlOt  OF 

THE  SPECTRUM. 

IF  LOW  IS  .FALSE.  I NO  CONVERGENCE  AT  THE  LOWER  SICE. 
-LOGICAL.  IF  Up  IS  .TRUE.  « CONVERGENCE  AT  THE  UPP£P  SICE  OF 

THE  SPEC  T RUN , 

IF  UP  IS  .FALSE.  * NC  CONVERGENCE  AT  THE  UPPER  SIOE. 
-LOGICAL.  IF  OIL  IS  .TRUE.  I EVSCAN  HAS  DcTERf'INEO  NDIV  fcIG£KVAlU£ 
INTERVALS,  WHICH  ARE  BELIEVED  TJ  REPRESENT  THE  NCIV 
S IALLEST  EIGENVALUES  OF  A (J.E.  ALL  NOIV  SMALLEST 
EIGENVALUES  OF  A HAVE  IEEN  DISCOVERED). 

CAUTION!  FROM  THE  NATURE  OF  THE  l A NC 70S - P R CC E SS  IT  IS 
CLEAR,  THAT  THIS  PARAMETER  HAY  YIELO  WRONG 
INFORNAT ION. 

-INTEGER.  SEE  OIV. 

-DIMENSION  TK(MT,2).  THE  °AIRS  ( TK  ( 1 , i ) , 7K.(  1 , 2 ))  , 

(TK (2, 1) ,TK (2,2  >)  , ...  , tTK(NtV.l) »TK(NEV«2)> 

REPRESENT  THE  NEV  INTERVALS,  WHICH  CONTAIN  EIGENVALUES 
OF  A. 

-INTEGER.  ERROR  INJICAT ION. 

THE  VALUE  OF  IE  R R IS  SET  EQUAL  TO  AN  CRRCR  COMPLETION  CODE, 
THE  NORMAL  COUPLET  ION  CODE  IS  0. 

IF  . 1 0 rl  _ THAN  3 J ITERATIONS  ARE  REQUIRED  TO  OETER  FINE 
AN  EIGENVALUE  (USING  THE  SUTROUTINE  IITQtl)  OF  T He 
SYMMETRIC  TRIOIAGONAL  MATRIX  T,  IERR  IS  SET  EQUAL  TO 
THE  I N DE  > OF  THE  EIGENVALUE  FOR  WHICH  THE  FAILURE  CCCUPS. 
-REAL.  A rIE  ASUR2  FOR  THE  RELATIVE  ACCURACY  IN  THE  CONFUTED 
EIGENVALUE  INTERVALS,  AS  DELIVERED  IN  TK. 


C INTERNALLY  CALLc.0  SUHPROGRA-S 

c .................. 

C IMTQL1  (A  SUBROUTINE  FROM  THE  E ISPA CK- aA CK AGE  1 

C MONO  IS 

C SCAN, IP 

C COUNT 

C EPSSP.N 

C (ALL  THESE  ROUTINES  ARE  INCLUDED  IN  THIS  DECK) 

C MLTPLT  (71520) 


C REMARKS 
C 

C I> 


EVSCAN  JAY  DE  CALLED  AFTER  A CALL  OF  LSVLJN  (70517)  OP 
GENLAN  (70513). 

SOMETIMES  EVSCAN  DELIVERS  AN  INTERVAL,  WHICH  CONTAINS  NO 
EIGENVALUE  OF  A.  HOWEVER  IN  SUCH  A CASE,  THERE  IS  AN  EIGENVALUE 
IN  Ttu  NE  1GHOOURHOOO  OF  THAT  INTERVAL.  IN  THESE  SITUATIONS  IT 
IS  CON  ION  THAT  THE  PROCESS  YIELOS  ALSO  Tilt  INTERVAL  IN  WHICH 
THE  RttSPECTIV.:  EIGENVALUE  IS  SITUATEO.  IN  ORDER  TO  IDENTIFY  BOTH 
INTERVALS  A 5 REPRESENTING  THE  SAME  EIGENVALUE,  IT  IS  ADVISED 
TO  USE  THE  SUBROUTINE  FLTPLT  (70520)  WITH  A SLIGHTLY  LARGER 
VALUE  OF  THE  OUTFUTVALUE  EPS  (10*CPS,  S A Y I t 
E PS  = 10 *EFS 

CALL  NTLFIT ( MT , NE V , EPS , TK I 

FROM  THE  NATURE  CF  TH  L ANC ZOS- PROCE SS . IT  FOLLOWS  THAT 
POSSIBLY  SOIL  EIGENVALUES  ARE  NOT  DISCOVERED  AT  ALL. 

IF  THE  ORIGINAL  MATRIX  A HAS  A MULTIPLE!,  THAN  THIS  MULTIPLICITY 
IS  IGNORED. 

EVSCAN  IS  WRITT-N  DY  J.M.  VAN  KATS  ANO  H . A . VAN  OCR  VOKST 
( ACCU)  UTRECHT. 

^ISP^ISBEStQOALmm^^ 

CUT  x 1 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 


C •***.»* 

C EXA  1PLE 

C •»*•»*♦ 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


JTW11  COPY  FURNISHED  10  DDC 

OF  USE 


IN  THE  FOLLOWING  flocks  t 16  L ANC20G- 5 TE PS  ARE  F ER^OR.-EC  CN  A 21TH 
ORDER  MATRIX  (WHICH  IS  SHOWN  AS  THE  FATRIX  W21  + OF  WIL<INSON>. 
EIGENVALUES  OF  THE  MATRIX  W21*  ARE  COM°UrEO  DY  EVSCAN,  USING 
THE  RESULTS  OF  LSV  LAN. 

PROGRAM  LNCZOS(CUTPUT) 

D LIENS  ION  QJ(2*),Q..(<l),U(2l)  , ALPHA ( 16) , PET A (16J 
LOGICAL  RSSTRT 
EXTERNAL  W21PLS.SAVEQ 
DIMENSION  TK  (16,2) ,TK1(16,2) 

LOGICAL  LCW,U°,OIV 
N = 21 

IFIRST=i 
U = 16 

RESTRT  =. FALSE. 

NSEEO=il J777 
CALL  RANSET ( NSEEO) 

DO  10  1=1, N 

10  Rid)  = RANF  (SEED) 

CALL  L3VIAN(N,IFIR3T,N,RESTRT,R1, W21PLS,SAVER,U,00 
♦ .ALPHA, BETA) 


flT  = 16 

CALL  EVSCAM.l,N,-T,ALoHAt3ETA,NEV,L0W,UP,DIV, 

♦ NCIV.TK.TKl, IERR.EPS) 

IF(IERR.NE.O)  STOP  "ERROR  IN  E IGE NVALUEC CM  PUT AT  ION" 

IF(LOW)  PRINT  1000 
IF (UP)  PRINT  1010 
IF  (O I V ) PRINT  1020.NDIV 
PRINT  1020, NEV 

C IF  NEV  EQUALS  0 NC  EIGENVALUE  INTERVALS  HAVE  BEEN  OETECTEO, 
C LSVLAN  SHOULD  eE  RESTARTEC  F OR  FURTHER  DETECTION,  WITH  THE 
C DIMENSIONS  OF  ALPHA,  BETA,  TR  AND  TK1  PROPERLY  AOJUSTSO. 


IF  (NEV.E 

C.O)  GOTO  70 

OO  60  1=1 

»NE  V 

60 

PRINT  10 ,0 , I ,TX (I , i) ,TK(I,2) 

70 

CONTINUE 

1000 

FORMAT (* 

CONVERGENCE  AT  LOWER 

SIDE  M 

1010 

FORMAT (* 

CONVERGENCE  AT  UPPER 

SIDE  *,/) 

1020 

FORMAT ( ♦ 

* , I 3 , * INTERVALS  AT 

LOWER  END  *) 

1GJ0 

FORMAT  (• 

*,//,*  *,I1,*  INTERVALS  ARE  FOUNC  *,/, 

• 

NR  LOWERBOUNO 

UPPEReCUND*,  /) 

1U0 

FORMAT  (* 

* , 1 3 , * *,2(E20.13,* 

*>  ) 

E NC 

* 

SUBROUTINE  SAVEC(Q.N) 

DIMENSION  Q (N) 

C NO  EIGENVECTORS  APE  COMPUTED,  SO  WE  DO  NOT  STORE 
C THE  ORTHOGONAL  TR A NSFOR 1 AT  ION- F ATR IX . 

RETURN 

ENO 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c ***.*.. 

C METHOO 

C 

c 

c 

c 

c 

c 

c 
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subroutine  hji  pl  3 ( v , ax,  rn 
0 IMENS  ION  X (N)  , fix  (N) 


THIS  PACK  IS  BEST  QPM  1 • ' ^ *vl' 
JTWJi  CcXi  T fUKMSilSl'  i*-'  " 


'U'Jjii 


C MATRIX  IS  TAKtN  FRO  I WILKINSON  (THE  ALGEBRAIC  £ IG  r N V A LU  £ PR  0 EL  c *■  ) 
C ANO  13  CALLED  W’l». 


C 

c 

c 

c 

c 

c 

c 

c 


10 

1 

0 


1 0 

BIO 
1 t 1 0 


0 

1 

10 


AX (1)  = 10»X (1) *X  (2» 

DO  10  1=2,23 

AX  (I)  = X (1-1)  ♦ IABS  (ii-IJ  »X  (I)  *X  (I  + U 
10  CONTINUE 

A X ( 2 1 1 = X ( 2C  J *13  *X (21) 

RETURN 
E NC 


THE  OUTPUT  OF  THIS  PSCGRAY  ISJ 

CONVERGENCE  A7  ICWER  SIDE 
CONVERGENCE  AT  UPPER  SITE 

2 INTERVALS  AT  LOWER  END 


6 INTERVALS  ARE  FOUNC 
NR  LOWERBOUNO  UaPERBOUNO 


1 Il2544ie:2il92*0  1 

2 «2533u53*7'JB7JT*J0 

3 . 7u03  )S22J9-,22:.»Jl 

4 . 3J33Bh1I223’6E>01 

5 . BJ10o736473 37E+01 

6 . 1074olS41d2B JE+t 2 


-.11254-,1522119E*-01 
•25TT058i70996E»CO 
. 70  33B5220<3-<24E*C1 
. 8J33B41l22  37bE»-Cl 
. B21j678o-<7277E*01 
.107431 94L329QE*02 


DISTURBANCE  OF  THE  INTERNAL  MONOTONY  OF  THE  EIGENVALUES  OF  SUCCESSIVE 
TRIOIAGONAL  MATRICES,  AS  0ESCRI3E9  INI 

VAN  KATS  J.N.,  VAN  OIR  VORST  H , A . , 

"AUTOMATICAL  fONITORING  OF  GENERALIZED  SYMMETRIC 
OR  SKEWSYf'iETRIC  LANCZOS  SCHEMES", 

1977,  A CCU,  UTRECHT,  T R 7 • 


SUJUD'JTINl  iVSCAN(H,  N,  M,  ALPHA,  ^TA,  U £V,  LOW,  11°,  CIV,  NDIV, 
i IK,  TKi,  IwRR,  E 3 S ) 

DI  I..NSION  ALPHK'O,  iUjai!)  , TK(,-T,2),  T K i ( :i  T , 2 ) 

L C » I C A L LCW,  U’,  01 V 

NiV  = 0 

LCW  = . F4L32. 

DP  = . FALSE, 

NO  IV  = 0 

io  r=2,  i 

J - .1  - I ♦ i 
T <(.)♦!,  2)  = 1 £ T 4 ( J ) 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
FROM  Cur  I FURNISHED  TO  DOC  -- 


.) j ::  i =i , i 

TK  < I, i>  - iL’HI ( I ) 
7N_  (1,1)  = AlcH4  ( I) 

:<_(:,£)  = r < ( i , ’ ) 

C u ■ ■i » I S U .. 


IF  (liRR.NE.J)  RETURN 

t<(  i, = : . 


£****»*»•**»*»»*»**•***»•*****»,**»»*¥■»** 

CALL  Ir.rJUCil,  TK,  T K ( 1 , 2 ) . IIRR) 

£»*»¥»**¥*»»¥*»*****¥•«•***<►*»*¥*****♦¥*»* 


£ FS  = FLOAT  (M  *AMAKi(CPS  ,2.  **  (-47)  ) 

C IF  CL".  WANTS  A LARGER  VALUE  CF  EPS,  THIS  VALUE  SHOULO  OE  INSERTED  HERE 
00  J j 1=*  , li 

T<(I, 2)  = TK( 1,1) 

TKi (1,2)  = TKi ( I , i ) 

3 J continue 

TKi  (.1,2)  = T K .,  { !,i) 

K = ■<1 
Kx  = l 

»¥*»»»»»»*»*»,»»»»¥*¥»»»»¥*»»*»»**♦»*»*»»*******»¥****¥,  ¥**,**»¥** 

CALL  CCA  N rP ( MT , K,  <1,  TK,  TKi,  EPS) 

C*  * * ¥**»***»*»*,»***»»»**¥*»**»**»,**»»¥*,,»*¥♦♦¥,¥*¥**,»¥¥¥¥,  *********, 

,j  CONTINUE 

E =S  OLD  = £ 33 

Q****»¥»*»»*,  *♦».»,  **»»»***¥*,***»¥*»***»,,*¥¥,*»******,*,»**, 

CALL  C 0. 1 1 N T ( .I , TK,  TKi,  K,  Ki,  E^S,  LON,  U°,  NEV,  DIV,  NCIV) 

IF  (EPS.GT.EPSOLO)  GO  TC  *0 
IF  (..NOT. LOW)  OIV  = .FALSE. 
return 


J 


FROM  CU1  \ FUKNXSHfiD  10  LOO  


S L»  ‘ i -J  C ‘ i T I N i.  S C A I i ? ( 1 , K , Ki( 

oi  i_w.»ijin  r<(.  tki  (/.ei 


TK,  T<;,  EPS) 


c with  the  val’J-  or  two  rcjs  (each  afarti  are  scanned  for  -ulticlet 

C T(IS  AY  A E L I V Z R 1 ’JEW  VALUE  Jc  Ep3  AND  THE  PROCESS  IS  RE  pE  3 T £ 0 AS 


c 

L 

I J G 

AS  EpS 

ll.  . j i 

c 

T 

K S 

RE  -St J L T 

S I J TWO  lULT 

: PLE  T 

-FREc 

‘"•O^  0 « 

c 

1 

I S 

THE 

JU 

HER  CF  ROWS 

IN  TK 

AND 

tk:# 

c 

K 

I S 

TH.. 

■JU 

HE  R OF  I N T E R 

VALS 

IN  TK 

• 

c 

< 

TrH  N 

1 4EER  OF  INTE 

RVALS 

IN  ~ 

Kl  . 

_3  : o ir i ju  : 

KOLD  = K 
K.OLQ  = K_ 

EpKOLD  = Zp3 

CALL  ILTPLTC.I,  K,  £?SOLD,  TKI 

call  iirnrr  ■,  <1,  eps,  r<i> 

EPS  - A I Ax;.  ( Z -3  )LC  , EPS) 

IE  ( (K..JE.KOLC)  .OR.  (Ki.NE.KiCLD)  ) GO  TO  10 

RET  )R  4 

END 

3UOJOJTINZ  C 0 • I'U  t I , TK  , ~Kl,  K,  Kl,  EPS  , LOWER,  upder,  jcee. 

1 OIV,  J 3 I 0 ) 

DIMENSION  TKt  ,2),  TK.-.l  r,2> 

LJfSICAL  L CNE  R , NPPER,  CISJCT 
LOGICAL  OIV 

C CO  I c A VISION  OF  TWO  ROWS  OF  INTERVALS  GIVEN  IN  TK  AND  TKI. 

C IF  AN  INTERVAL  IN  Oil  ROW  IS  CLOSE  TO  AN  INTERVAL  IN  THE  OTHER  ROW. 

C THAN  THE  jpAN  OF  4 0 7 H INTERVALS  IS  RECCRCED  IN  TK. 

C SO  THE  VALUES  OF  7K  HAVE  DEEN  0 VE  R WR  I T T E N . 

C LOWER  IS  .TRUE.  : CONVERGENCE  AT  THE  LONER  END  OF  THE  SFECTfiUN. 

C IJ 3 IS  .TRUE.  : CONVERGENCE  AT  THE  U pF£  R END  Cc  THE  SPECTRUM. 

C OIV  IS  .TRUE.  » A HOLE  IN  THE  S c E C T A L m IS  DETECTED. 

C J 3E  F IS  THE  Nl  IDE  R OF  THE  CC-'PUTED  EIGENVALUE  INTERVALS. 

C J II  V GIVES  THE  .JU  HER  Op  EIGENVALUE  INTERVALS  AT  THE  LOWER  ENC. 

OIV  = .FALSE. 

IOIS  = { 

LOWER  - .FALSE. 

UPPER  = .FALSE. 

FOLO  = A UNldKld.D.TKtl.i)) 

EOLJ  - (l.-SIGNC  j.l.FOLOM  *FOLC 
FOLD  = EOLO 
J = 1 
J-  = 1 
JD-EF  = 0 
-0  CONTINUE 

CALL  EPSSFNt  TK  CJ  , 1)  , TK  ( J » 2 ) . TKUJl.il,  TKKJ1.2),  EPS,  E,  F, 

_ DISJCT) 

IF  (CISJCT)  GC  TO  53 
I D I S = j 

IF  (J.EQ.i  .OR.  Ji.EQ.il  LOWER  = .TRUE. 

IF  (J.EQ.K  .OR.  J 1 . E Q . K 1 I UPPER  = .TRUE. 

IF  (F.LE.FOLO)  GO  TO  EJ 

IF  (E.GT.tOLD)  J D E f = JCZC  + 1 

T K ( J DE  c , 1 ) = E 

£ OLD  - E 

TK(  IDEF, 2)  = c 

FOLO  = F 

23  IF  ( r< (j, 2) . L£. TK_ ( Ji , 21)  GO  TC  -*0 
JJ  CONTINUE 

J = J f 1 

IF  CJ.LE.K)  GG  TO  iu 


A 


ISKfW 


- A3 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
EROJii  COPY  FURBISHED  TO LUC 


I = J - 1 
J*  = Jl  *■  1 
I c r.<i)  go  to  a: 

,!•  CONTINUE 

J » = I i ♦ J. 

IF  ( J.iLi.Ki)  jO  TO  1 u 

Ji  = ji  - * 
j = j ♦ i 

I T ( J . G T . K ) GO  TO  .*} 0 
>j  C 0 J 7 Ii  U E 

C IF  SUCCE  S 3 IVE  l Y 3 0I5JUNCT  INT£qVALS  ARE  'ET  , A HOLE  IN  THE  SPECTRUM 
C II  SIJF303EO. 

ICI3  = 1313  1 

IF  (.'JOT. CIV  . A ' 1 0 . (IOIS.EQ.c)  ) GO  TO  60 
GO  TO  ?\j 
6j  JOI  V = ,13 : F 
3j.V  - • T R L ^ • 

70  IF  (T<i(Jj.,i)  .GT  «TK  (J  ii)  ) GO  TC  3G 

GO  TO  *t 0 
30  CONTINUE 

I F (.NOT. (LOWER  , AND . (.NOT.OIV)))  RETURN 
0 1 V = .TRIE. 

JOI/  - JOcF 

R-T'JRN 

ECO 

SU3  ROUT  I Nc  E=  SS  = MA,  .3,  C,  0,  EPS,  E,  F,  OISJCT) 

LOGICAL  OISJCT 

C IF  TWO  INTERVALS  IA,3]  AUO  I C , 3 1 ARE  RE l AT I VEL  Y CLOSE  WITH  RE  S°ECT  TC 
C EPS,  ( E , F ] GIVES  THE  SPAN  OF  GOTH,  AND  DI3JCT  IS  SET  TO  .FALSE. 

C IF  THEY  ARE  NCT  CLOSE,  OISJCT  IS  SET  TO  .TRUE.. 

■DIS.ICT  = .TRUE. 

e °ss  = :. 

IF  (A.LE.C)  GO  TO  10 
IF  (O.GE.A)  GC  TO  20 

IF  ( A 3S(  3-A  > / (1.  ► AOSCO)  ) .LE.E5S)  GO  TO  20 
RET  PEN 

lO  IF  (C.LE.PJ  GC  TO  23 

IF  (A3StC-n>/(l.+A93(C> ) .LE.EPS)  GO  TO  23 
RETURN 
20  CONTINUE 

OISJCT  = .FALSE. 

E = A Hill  ( A,  C) 

F = A I A X 1 (13,  C) 

I PSS  = (F-E)/ti.  ♦ A>OS<E)  ) 

IF  (EP3G.GT.EFS)  EPS  = EF3S 
RET  J <N 
E NO 

SUG  ROUTINE  IOiNOIS  (TK,  TKl,  K,  E p 3 3 
01  lENSION  TK  ( K ) , TK_  ( _) 

REAL  LON 

C OIJTURGANCE  OF  ,'lONQTONY  OF  TUO  SUCCESSIVE  RCWS  OF  EIGENVALUES  GIVEN  IN 
C A NO  TKL. 

C TKl  HAS  O J£  EIGl  (VALUE  .lO'RE  THAN  TK  (K  + l AND  K RES^.t. 

IJP  = T K ( 1 > 

VAL  = TKl  ( 1 3 

E PS  = j . 'J 

IF  (VAL.LE.UP)  GO  TO  iO 

EPS  = AGS (VAL-U^) / (AOS (VAL) ♦«. > 

*0  CONTINUE 


AD-A061  729 


UNCLASSIFIED 


ACADEMIC  COMPUTER  CENTER  UTRECHT  (NETHERLANDS)  F/G 

FAST  ITERATIVE  METHODS  FOR  LAR6E  LINEAR  SYSTEMS. (U) 

AUG  78  H A VORST  OA-ERO-75-G-0084 


F/G  12/1 


2°f2 

£861729 


TO  i J I = 2.K  THIS  Pi 

L T I - f fMOU  Ct 

VAL  - TKi ( I ) 
l=>  = TKO 

if  (val.g:.ltn>  go  rc  2j 
_ = = AES f VAL-LCW) / (AES  (VAL > ♦! . > 

(I  = .GT.  ;oS>  EPS  = EP 
If  <VA L.LE.NJ>  GO  TO  ?J 
_P  = A63(UF-VAL)/(A-1S  (VAD+l.) 

IF  (ZP.GT.  i=>Z)  ZPS  = £3 
CONT  I VUE 
VmL  = TKiK  + i) 

IF  (VAL.  Gc.U3)  GO  TO  j 

£3  = «1S  (VAl-LPi/tA  IS  (VAL)*!.) 

IF  (E°.GT,ZPS)  lc5  = £ P 
CONTINUE 

IS  THE  lAXI.lU.t  OF  THE  RELATIVE  GIST  A NCE 
. MONOTONY  CRITIRIUh. 

R IT  JRN 
£ MO 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 

f»U*l  CUi  i EUKMSHED  TO  DOC  ^ 


OF  OISTtREANCE  CF 


3lJ‘B  ROUTINE  I it  TO  Li  (N  » 0,  E*  IERR) 

INTEGER  I,  J,  L,  il , N,  II,  fihL,  IERR 
REAL  D < N>  , Z (N) 

REAL  0,  C,  F,  G,  =>,  R , S,  , I A CH  c P 
C RIAL  SO  IT ,A OS, SIGN 

C THIS  SUBROUTINE  IS  A TRANSLATION  OF  THE  ALGOL  PROCZOURE  IMCLi, 
C NUN.  iiATH.  l2,  077-3  33  (_96S)  OY  MRTIN  A NO  WILKINSON, 

C AS  hODIFIEO  IN  UU;i.  1ATH.  15,  A50(1B7G)  OY  OUORULLE. 

C HANDBOOK  FOR  AUTO.  C3>°.,  V C l . I I- L INE A R ALGE3RA,  2 A 1 - 24 8 ( 1 9 71 ) , 

C THIS  SUBROUTINE  FIDOS  THE  EIGENVALUES  CF  A SY, METRIC 
C T RI DIAGONAL  MATRIX  OY  THE  INFLICIT  OL  METHOD. 

C OU  INPUT- 

C N IS  THE  ORDER  OF  THE  ,;ATFIX, 

C C CONTAINS  THE  DIAGONAL  ELEMENTS  OF  THE  INPUT  MATRIX, 

C E CONTAINS  THE  SU3C IAGONAl  ELE  TENTS  CF  THE  INPUT  MATRIX 

C IN  ITS  LAST  N-,.  POSITICNS.  E(i)  IS  AROITRARY, 

C ON  OUT  PUT  - 

C C CONTAINS  THE  EIGENVALUES  IN  ASCENDING  ORDER.  IF  AN 
C ERROR  EXIT  IS  RACE,  THE  EIGENVALUES  ARE  CORRECT  AND 

C ORDERED  FOR  INDICES  1,2,... IE RR-1 , OUT  MAY  NOT  BE 

C THE  SMALLEST  EIGENVALUES, 

C £ HAS  OEEN  DESTROYED, 

C IE  RR  IS  SET  TO 

C ZERO  FZi  NORMAL  RETURN, 

C J IF  THE  J-TH  EIGENVALUE  HAS  NCT  BEEN 

C DETERMINED  AFTER  3D  ITERATIONS. 

C QUESTIONS  AND  CO  i.'ZNTS  SHOULC  BE  DIRECTED  TO  B.  S.  GARBOW, 

C APPLIED  MATHEMATICS  DIVISION,  ARGONNE  NATIONAL  LABORATORY 

C 

C **********  iACUEP  IS  a MACHINE  DEPENOENT  DARAM£TER  SPECIFYING 
C THE  RELATIVE  PRECISION  C F FLOATING  °OINT  ARITHMETIC. 

C ********** 

IACHEP  = 2.**<-*7> 

I ERR  = D 

IF  (N.  £0.  II  GO  TO  l4i) 

DO  1=2,  J 

E(I-i)  = E(I) 

1C  CONTINUE 


J 


c 


c 


c 

c 


3 j 

j 


60 


7 J 


10 


li> 

LJU 

lij 
12  J 
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THIS  PACK  *3  BEST  QUALITY  PRACTICABLE 

£ Cl)  - J.3  FROM  COl  Y KUKMSIUiU  J\)  Ui»C  

1C  _J  L = A,.'I 
J - , 

******  LOC<  FO  » ALL  SliE-OIAGCNAL  cLcMEAT  ********** 

10  5 J f*  = L • 

IF  «.  I.  :q.:i»  GO  TO  BJ 

IF  (AHS(£(  1>)  ,LC..*  flCHFP*  (AOS  (C  ( 1)  > +A1S(D(«*1>  > >>  GO  TO  4C 

GONT  INLC 
? = 1 ( L ) 

IF  ( I.  ; L)  GO  TO  Id 
IF  ( ).  - ).  3J)  GO  TO  12d 
0 = J ♦ 1 

******  FOB  1 SHIFT  ********** 

g = (jiLu)-3)/(:.j*e(in 
•i  = SORT (G*G*i.G) 

G = 0(.')  - " * E (L)/(G*SIGM(*»G)> 

5 - * . u 

L*  ~ x • 0 

p = 0 • u 
. f :i  L = r - L 

'******  F O h 1 = --.  STEP  -1  UNTIL  L 00  — ********** 

10  7 J 1 1 = x » I )L 
I = .1  - II 
F = S V (I) 

1 = 0*~  ( I> 

IF  IA3S(F)  .LT.A-3SCG)  ) GO  TO  50 
C = G/F 

R = S.KT  (C*C*i.  U) 

z ( I *1)  = r*S 

S = 1.0 /R 

c = c*s 

GO  TO  oO 
S = F/G 

R = S )RT ( 3*3*1. J) 

£ ( 1*1)  = G*R 

c = L.d/R 

3 = 3*  C 

G = 0(1*1)  - o 
R = (0 ( I ) ~G) *3  ♦ c .0  *C* 1 

P = S*R 

0(1*1)  = G ♦ ° 

g = c*r  - a 
COOTlNLt 
J(L)  = D(L)  - P 
£ ( L ) = G 

£ ( ! ) = J . 0 
GO  TO  20 

******  ORCER  EIGENVALUES  ********** 

IF  (L.cl.l)  GO  TO  jlJO 

******  FOB  I = L STEP  -1  UMIL  2 00  --  ********** 

00  10  1 1 =2  *L 

I = L ♦ 2 - II 
IF  (P.Gc.O(I-l) ) GC  TO  llv 
0(1)  = D(I-i) 

CONTINUE 

1 = i 

D ( I ) = p 
CON T INUE 


ir 
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THIS  V'AQE  . wuuC 

n*« 


j ro  t .j 

c **••****••  j_j  0\ji  — uo  l'inv  -?r.£NC-r  ir  am 

C £ ICw  IV  ^ L ’ J £ 4 r T £ S cl  I T " F A T ICNS  * * ‘ 

i ij  : n*  *.  = l 

l .j  .R.T  K'  U 

C * L A '»  7 CA-’l)  OF  i t T : L 1 •******•< 

c tat:  j./;j/76 

C *****••_  40  OF  J £(<*•****•***  I « T f;  L I 


HEADING  fuSCl) 


SUBROLT  III.  ILIPL  r C.'I.MNT.c  F S,  A I 
01. It  NS  ION  A ( ,'t  t 7 ) 


PUR  FOSE 


Ul  KEPLACl,  IN  A <OH  or  IMlRVALS,  THOSE  INTERVALS  TMAT  ARE  RELATIVELY 

ct  o ;t:  iy  m ir  span. 


INPUT  PARAMETERS 


-INI  !(.£’{.  lilt  Nil.' HE  R OF  ROWS  IN  THE  All  Til  A L OE  CL  AR  AT  TON  CF  THE 
ARRAY  A. 

-INTeGlR.  T Hi  Nil  FOE  R OF  INTcJVALS  GIV  N IN  ARRAY  A. 

(ALSO  OUT  HIT  PAN  AM  til 

-REAL.  IF  Tilt  OISTAUC.  Of  TWC  SUCCESSIVE  INTERVALS  IS  RflATIVIlY 
LESS  THAN  tPS,  THEY  ARi  Re  PL ACC 0 BY  THEIR  SPAN. 

IALSO  OUT  FU  T PARA  tv  T si 

-OIMENSION  AIi.JT.  IH  pairs  IAlMI,A(l,:il,  C AI?,1»  .AIC.2TT  . 

I AININI  , l»  , A _ p R • S N T TH  M NT  INTERVALS. 

THE  FOLLOWING  CCNOITICNS  SHOULD  HCLOl 
i . A ( I , J ) . f.  £ . A ( I . 1 1 
?.  A (I  ♦!  ,11  , GE  . A ( I ,;i 

(ALSO  OUIFIIT  PARAMETER! 


OUT  FUT  PARA.lt.  T.  RS 


-INTEGER.  IMF  NUFO.R  OF  OISJUNCT  INTERVALS  IN  TMF  FINAL  ROW  A. 
(ALSO  IMUT  ‘•ARAnETERI 

-REAL.  THE  NAYINL  I OF  THE  RUATIV  WIOTHS  OF  THE  INTERVAL?  IN 
THE  FINAL  ROW  A,  0. F I N 0 «V| 

NAM(A(I,E>-a<r,l)>/AAS(4<r,l>Rl.0>> 

IF  THIS  VALUe  IS  SlALLER  THAN  THE  I N FUT  VALU  OF  E BS  THFN 
EPS  IS  NCI  CHANCE 0. 

(ALSO  I N I LIT  PA  ‘A  1 1.  T E N I 

-OITIeNSION  AIn.’I.  CONTAINS  TMl  R SULTING  NINT  INTERVALS. 

I ALSO  INFUT  PARA.iETF 


INTERNALLY  CALLt n SUOPROl.  TANS 


Rt H ARNS 


Ml  T PL  T tAY  BE  U Sc  0 AFT.R  EVSCAN  (.'JSlBI  IN  OR  OF  R TO  R MOVE 
POSSIBLY  S>»’IRIOUS  R SUITS. 

FOR  0.  TA  ILS  SEE  I 

VAN  NATS  J.N.,  VAN  0* R VORST  M.A., 

"Autonat  icai  mon  i t or  I Ni  or  general  i.’fO  sympftnic 


r 


os  SKEWSYN.' ETRIC  LANC70S  SCMEMcS", 

1177,  ACCU , UTRECHT,  TP  7, 

Ml  T PL  T MAS  TEEN  WRITTEN  HT  J..(.  VA  N KATS  l M.A.  VAN  Of!*  VORST 
( ACCUI  - UTRECHT 


txA.iPit  or  USl 


PRO CRAM  Ml TT  £S (OUT  PUT  1 
DIMENSION  A(^,2t 
rt3  *♦ 

NI  NT  3 *♦ 

EPS3i.E-2 
A(l,l>  «i>.  3 
All, 2)  3C. 3 0 5 
A(2,l> *0. 51 
A (2 ,2)  3 j. 5 i 
A ( 3 , 1 ) =0.S20S 
A(J,2»  =0.53 
A ( <«  ( 1 ) =0.  71 
A(N,2>  =0.77 

CALL  ML TPl T (l,NINT,t PS,A) 

PRINT  1E30  .NINT.EPS 

for.;at(*  •,•  nu  'per  of  final  intervals 

• • EPS  *,E1.3,//»  . 

PRINT  1„1J 

DO  11)  I31,N  INT 

PRINT  l02u.A(I,l),A(I,2l 

FOR, 'AT  ( • »,•  IOWERHOUNO  IIPPSROOUNO  •» 

FORMAT ( * *,2(£il,bt2X)l 

END 


*, 12 ,/« 


THE  OUTPUT  OF  THIS  PROC.RAM  IS« 


NUrINCR  or  FINAL  INTERVALS 
EPS  .*J2E-Jl 


LOME  RflOUNO  UPF  t RDOUNC 

• 3u Ju JO  .♦  C J •3jiJOO-.*03 

.51t)«00i*0  J .53JJOEE*00 

• 76J0D0E  *00  . 77JJ»it't*03 


S(Jf3  ROUT  I Nt  ■ IL  T r*L  T ( ,1  , NIM,  EPS  , A> 

DIMENSION  A(.M,2> 

C THE  PART  A (NIM,  2)  CONTAINS  THF!  MUL  T I^LE  T-I NTE  R VALS  . 

CPS  I = CPS 

I r 1 

C J l.  IV.  S THc.  NLHiKK  OF  THE  CURRENT  MJLT  JPLCT . 

OO  E0  I -E  » N IN  T 

IF  (AMS ( A ( I , „) -A (J , 2 > )/ < 1 . ^AHS ( A ( I, 1) ) ) .GT .EPS)  GO  TC  10 
A(.I,E)  = A (I  ,E) 

C THt  UPP.rWMOUNE  i)F  TMt  J-TM  .‘ULTIPLfT  IS  UPDATED. 

EPS, IUL  = AOS  ( A (J  ,E  ) - I ( J,  l>  ) 7(  1.  ♦ APS  ( A ( .1,1)  ) ) 

IF  C PS-U.IL  .G  T . PSM)  £ f S 1 = EPS-Ill 
GO  TO  2U 
iJ  - J ♦ i 

A (.1,1)  = A (1,1) 

A ( J , .?  ) = A (I  ,2) 

20  CONTINUE  .^aOTYCNVA** 
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m 

number  of 
eigen- 
values 

BP 

€ 

H 

— 

number  of 
eigen- 
values 

— 

m 

H 

k 

Mm 

H 

21 

40 

152. 

7.2 

22 

40 

132. 

6.0 

3 

0 

i 

0.3 

23 

40 

276. 

12.0 

4 

24  >) 

5.1 

1.3 

24 

40 

188. 

5 

32  >) 

21. 
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25 

40 

156. 

6 

35  ’) 

20. 
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26 

40 

159. 

6,1 
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40 

13. 
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27 

40 

287. 

10.6 

a 

40 

44. 

5.5 

28 

40 

275. 

9.8 

q 

40 

33. 

3.7 

29 

41  ;) 

224. 

7.7 

10 

40 

6?. 

6.2 

30 

41  ') 

114. 

3.8 

1 1 

40 

*>6. 

5. 1 

31 

40 

262. 

8.5 

12 

40 

130. 

1 1 .0 

32 

40 

282. 

8.8 

13 

40 

93. 

7.2 

33 

40 

281. 

8.5 

14 

40 

72. 

5.1 

34 

40 

192. 

5.6 

15 

40 

117. 

7.8 

35 

40 

359. 

10.2 

16 

40 

lib. 

7.3 

36 

40 

167. 

4.6 

17 

40 

61. 

3.6 

37 

42  ?) 

388. 

10.5 

18 

40 

183. 

10.2 

38 

40 

346. 

9.1 

19 

40 

205. 

10.8 

39 

40 

228. 

5.8 

20 

40 

288. 

14.4 

40 

40 

215. 

5.4 

Table  I 

Results  of  the  monitoring  process  for  the  har  prohlem. 


')  for  k«4,5,6,  all  eigenvalues  were  detected  at  the  upper  end  of  the 
spectrum. 

■)  if  the  recommended  extra  scan  is  performed,  40  eigenvalues  are 
detected. 


J 


SI 


-1.1254415221  19 
0.2538058170974 
0.9475343675298 
1.789321352695 
2.130209219363 
2.961058884186 
3.043099292579 
3.996048201383 
4.004354023440 
4.999782477742 
5.000244425001 


6.000217522256 
6.000234031583 
7.003951798616 
7.003952209528 
8.038941115813 
8.038941 122828 
9.210678647304 
9.210678647360 
10.74619418290339 
10.74619418290332 


Table  II:  Eigenvalue  of 


+ 10.74619418290 
+ 9.210678647331 

+ 8.038941119304 

+ 7.003952002664 

+ 6.000225680184 

+ 5.000008158672 

+ 4.000000205070 

+ 3.000000003808 

+ 2.000000000054 

+ 1.000000000000 
0 

Table  III:  Eigenvalues  of  j 


in 

number  of 
eigenvalues 

eps 

eps 

21* 2-47 

13 

0 

. I5E-I2 

1 

16 

6 

. 1 3E- 1 1 

8.9 

32 

20 

.83E-12 

5.5 

150 

19 

. 19E-10 

129. 

300 

19 

. 16E-10 

109. 

Table 

IVa:  Results  of 
the  Lanczos 

EVSCAN  after  m 
-scheme  applied 

iterations  of 
to  W*  . 

ra 

number  of 
eigenvalues 

ops 

ops 

2 l'*2~47 

13 

0 

. 15E-12 

1 

16 

0 

. 15E-12 

1 

32 

21 

. 1 5E- 1 2 

20.9 

150 

21 

.31E-11 

93. 

300 

21 

1 

o 

196. 

Table  IVb:  Results  of  EVSCAN  after  m iterations  of 
the  l.anczos-scheme  applied  to  W~  . 

I 


No 

lowerbound 

upperbound 

upperbound- 

lowcrbound 

1 

1 1254415221 19K+01 

-.  1 1254415221 19E+01 

. 135E-12 

2 

. 2538058 1 70973E+00 

.2538058 170981 E+00 

.815E-12 

3 

. 7003952209096E+0 1 

. 7003952209 I00E+0I 

.38  IE- 1 1 

4 

.8038941 1 1 9 703E+0 1 

.8038941 1 1 9703E+0 1 

. 1 14E-12 

5 

.92  10678647329E+0 1 

.92 1 0678647329E+0 1 

.568E-13 

6 

. 10746194 1 8290E+02 

. 1074619418290E+02 

0. 

Table  Va:  m-16 

No 

lowerbound 

upperbound 

upperbound- 

lowerbound 

1 

-. 1 12544 1 522 1 1 9E+0 1 

-. 1 1254415221 I8E+0I 

. 104E-1 1 

2 

. 2538058 1 70977E+00 

.25380581 709 80E+00 

. 293E- 1 2 

3 

.9475343675301 E+00 

. 9475343675303E+00 

. 263E- 1 2 

4 

. 178932 1352695E+01 

. 1789321352695E+01 

. 284E- 1 3 

5 

. 2 1 302092 19363E+0 1 

. 213020921936 3E+01 

.71  IE- 13 

6 

.296 1058884 185E+01 

.296 1058884 1 85E+0 1 

. 156E-12 

7 

. 3043099292578E+0 1 

. 3043099292578E+0 1 

.568E-13 

8 

. 399604820 1383E+01 

. 399604820 1 383E+0 1 

. 1 14E-I2 

9 

. 400435402344 1 E+0 1 

. 4004 3540 2 344 1 E+0 1 

0. 

10 

.499<?782477742E+01 

.4999782477742E+01 

0. 

1 1 

. 500024442500 1 E+0 1 

. 500024442500 1 E+0 1 

.853E-13 

12 

.60002 1 7522256E+0 1 

.60002 17522256E+01 

. 171E-12 

13 

.6000234031 583E+0 1 

.600023403 1584E+01 

.853E-13 

14 

. 700395 17986 I5E+0 1 

.7003951 7986 15E+01 

.256E-12 

15 

. 7003952209528E+0 1 

. 7003952209528E+0 I 

. 1 14E-12 

16 

.8038941 1 15813E+01 

.8038941 1 15813E+01 

.568E-13 

17 

.8038941 I 22828E+01 

.8038941 1 22828E+0 1 

. 1 14E-I2 

18 

.92 10678647304E+01 

.921 C678647304E+0 1 

0. 

19 

.92 10678647360E+0 1 

.92 1 0678647360E+0 1 

0. 

20 

. 107461941 8290E+02 

. 107461941 8290E+02 

. 125E-1 1 

Table  Vb 


m-32 


No 

lowerbound 

upperbound 

upperbound- 
1 owe r bound 

1 

1 1 2544  15221 18E+0I 

1 1254415221  1 1E+0I 

. 7 1 4E- 1 1 

'J 

.25380581 70894E+00 

.25  38058 17101 6E+-00 

. 122E-I0 

3 

. 947534  36753 1 2E*00 

. 947534 3675369E+00 

. 575E- 1 1 

4 

. 178932  l35269bF.t0l 

. 1 789 32I352700E+0I 

. 442E- 1 1 

5 

. 2 1 302092 1 9362E+0 1 

.21  302092 1 9364E+0 1 

. 205E- 1 1 

6 

.296 1058884 184E+0I 

.29610588841 86F.+-0 1 

. 189E-I 1 

7 

. 30  4 3099292573E+0I 

. 3043099292586E+01 

. 794E-1  1 

8 

. 399604820 1 379E+0 1 

. 399604820 1384E+0I 

.493E-1 1 

Q 

. 4004 35402 34 35E+01 

.4004 35402 3442E+01 

. 6 1 4E- 1 1 

10 

.49997824 7 7 739E+01 

.4999  78247774  1F+01 

. 350K-1  1 

i i 

.5000244424997E+-01 

. 500024442500 IE+01 

. 4 1 5E-  1 1 

12 

.6000 2 1 75  22  252E+0 1 

.60002  175222561:  HU 

.449E-1 1 

13 

.60002 3403I577E+0I 

.600023403  1 583K+-0  1 

.631E-1  1 

14 

.700395179861  1F.+01 

. 700395 1798763E+01 

. 152E-09 

IS 

. 7003952209506E+0 1 

. 7003952209528E+0 1 

.218E-I0 

16 

.8038941 1 1 5805E+0 1 

.8038941 1 15813E+01 

. 767E-I  1 

17 

.803894 1 1 2282 IK+01 

.80  38941 1 22828K+0 1 

. 682E- 1 1 

18 

.92106  78647298E+0 1 

. 92  lib'  7864 7 3b0EH)  1 

.619K-I0 

Id 

. 107461941 8289E+02 

. 107461941 8289E+02 

.483E-1  1 

Table  Vc : m=150 

No 

lowerbound 

upperbound 

vipperbound- 
1 over bound 

1 

-.  1 1254415221  19E+01 

-.  1 125441 522 106E+01 

. 128E-10 

J 

.25 380581 7090 1K+00 

.25380581 7 104 1E+00 

. 140K-10 

3 

.9475 34 367532 7E+00 

.9475 34 3675446E+00 

. 1 19E-10 

4 

. 1 789  32 1 352690E+0 1 

. 1 78932 1352702E+01 

. I26F.-10 

5 

.2 13020921 9 362E+01 

.2 13020921 9 372E+01 

. 958F- 1 1 

6 

. 2961058884 1 36E+0  1 

. 2961 058884 187E+01 

, 510E-1Q 

7 

. 3043099292577E+0 I 

.304 309929259 1 K+0 1 

. 1 39E-I0 

8 

. 399604820 1 372E+0 1 

.399604 820 1384E+01 

. 1 12E-10 

9 

.40043540234 34E+01 

.400435402344 IE+0 1 

. 685E-1 1 

10 

. 4999  7824777  37K+0 1 

.4999782477745F.+01 

. 867E-1 1 

1 1 

. 5000 2444 2499  3E+0 1 

. 500024442500 1K+0 1 

.801F.-1 1 

12 

.600021 7522246E+01 

.60002 1 7522256E+0 1 

. 10  IE- 10 

1 3 

.60002340 3 1570K+01 

.60002  3403 1 583E+0 1 

. I32E-10 

14 

. 700395 1798607E+01 

. 700  39  s 1 7986  1 4F.+0  1 

. 747E- 1 1 

15 

. 700  3952209507F.+0 1 

. 70039  52209  52  7F.+01 

.20  IE- 10 

16 

.8038941 1 15798E+01 

.8038941 11581 3E+0 I 

. 149E-10 

17 

.8038941 I22814E+01 

.8038941 1 22828E+0 1 

. 1 36 E- 10 

18 

.9710678647290F.+01 

.92  10678647  359F.+01 

. 697K- 10 

19 

. 10746194 18288F+02 

. 1074b  1941 8288F.+02 

. 750K-1 1 

Table  Vd:  in* TOO 


Eigenvalue  intervals  determined  by  KVSCAN  after  m l.anczos  iterations 

■f 

applied  to  W,. 


J 


< 
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No 

lowerbound  (a) 

upperbound  (b) 

b-a 

i 

I0746194I8290E+02 

-. 10746  19418290E+02 

. 248E- 1 1 

2 

-.92 106 7864 7 332 E+0 1 

-.92  1 06 7 864 7 3 32 E+0 1 

•227E-12 

3 

-.8038941 1 1 9305E+0 1 

-.8038941 1 19305E+0 1 

. 1 14E-12 

*4 

-. 7003952002662E+0 1 

-. 7003952002662E+0 I 

.256E-12 

5 

-.6000225680 1 84E+0 1 

-.6000225680 1 83E+0 1 

. 142E-12 

6 

-. 5000008 1 5867 1 E+0 1 

-.5000008 158671 E+0 1 

. 227E- 1 2 

7 

- . 4000000 20 50 70E+0 1 

-.4000000205070E+01 

.227E-12 

8 

-.300000000 330 7E+01 

- . 300000000 3806E+0 1 

. 384E- I 2 

9 

- . 2000000000054E+0 1 

-. 2000000000054E+0 1 

. 185E-12 

10 

-. IOOOOOOOOOOOOE+OI 

-.9999999999998E+00 

. 298E- 1 2 

1 1 

. 753 1604 1 83839E- 1 2 

. 10896222 1 9482E- 1 1 

. 336E-1 2 

12 

. 1000000000001 E+0  1 

. 1000000000001E+01 

.39  IE- 12 

13 

. 200000000005 3E+01 

. 200000000005 3E+01 

•284E-12 

14 

. 3000000003807E+0 1 

. 3000000003807E+0 1 

. 568E- 1 3 

15 

. 4000000205068E+0 1 

.4000000205068E+0 1 

. 568E- 1 3 

16 

. 5000008 15867 1E+0 1 

. 5000008 158671 E+0 1 

. 199E-12 

17 

.60002256801 8 1 E+0 1 

. 6000225680 1 82E+0 1 

. 1 7 1 E- 1 2 

18 

.7003952002661 E+0 1 

.7003952002661 E+0 1 

. 142E-12 

19 

.8038941 1 1 9 301 E+0  1 

.8038941 1 19301E+01 

0. 

20 

.92 10678647327E+0 1 

.92  1 06  786  4 7 32  7 E+0 1 

. 568E- 1 3 

21 

. 10746 194 18290E+02 

. 10746 1 94 I 8290E+02 

.210E-I I 

Table  Via: 

m=32 

No 

lowerbound  (a) 

upperbound  (b) 

b-a 

1 

-. 10746194 1 8290E+02 

-. 107461941 8289E+02 

. 144E-10 

2 

-.9210678647332E+01 

-.921067864731 9E+0 1 

. 131E-10 

3 

-.8038941 1 1 9 305E+0 1 

-.8038941 1 19291 E+0 1 

. 138E-10 

4 

-. 700395 200266 2E+0 1 

-. 700 39 5 200265 3E+0 1 

.901E-1 1 

5 

-.6000225680 1 84E+0 1 

-.6000225680 1 76E+0 1 

. 782E- 1 1 

6 

- . 5000008 158671 E+0 1 

-.5000008 1 58662E+0 1 

. 892E-1 1 

7 

-. 4000000 2050 70E+0 1 

-.4000000 205060E+0 I 

•955E-I I 

8 

-. 3000000003808E+0 1 

-. 3000000003799E+0 1 

.904E-1 1 

9 

- . 2000000000055E+0 1 

-. 2000000000046E+0 1 

.9 1 8E- 1 1 

10 

-. 1000000000002E+0 1 

-.9 9999 99 9999 2 4E+0 1 

.952E-1 1 

1 1 

-.  I337897763203E-I  1 

.4887256470544E-1 1 

.623E-11 

12 

.9999999999981 E+00 

. 1 000000000002E+0 1 

. 352E-1 1 

13 

. 2000000000050E+0 1 

. 2000000000054E+0 1 

. 355E-1 1 

14 

.300000000 380 3E+01 

. 3000000003806E+0 1 

•315E-1 1 

15 

.4000000205063E+0 1 

. 400000020506 7E+01 

. 406E-I I 

16 

. 5000008 158665E+01 

. 50000081 58667E+0 1 

. 222E-1 1 

17 

. 6000225680 1 76E+0 1 

.6000225680 1 82E+0 1 

, 585E-1 1 

18 

. 7003952002654F.+0  1 

. 7003952002660E+0 1 

. 560E-1 1 

19 

.80389411 1 9292E+0 1 

.8038941 1 19314E+01 

. 22  IE- 10 

20 

.92  106786473 1.5E+01 

.921067864731 9E+0 1 

. 38 1 E-l 1 

21 

. 107461941 8288E+02 

. 107461941 8289E+02 

.995E-I 1 

Table  VIb:  tn=150 


56 


No 

lowerbound  (a) 

upperbound  (b) 

b-a 

1 

10746 194 18290E+02 

-. 1074619418287E+02 

. 335E-10 

2 

-.92 1067864733 1E+01 

t.92  10678647302E+01 

. 290E- 10 

3 

-.80389411 19305E+01 

-.80389411 19280E+01 

.254E-10 

4 

-. 7003952002662E+0 1 

- . 7003952002643E+0 1 

. 196E-10 

5 

-.6000225680 1 84E+0 1 

-.6000225680158E+01 

. 26  IE- 10 

6 

-.5000008 15867 1E+0 1 

-.500000815865 1E+01 

. 206E-10 

7 

-. 4000000205070E+0 1 

- . 40000002049  3 1 E+0 1 

. 140E-09 

8 

-. 3000000003809E+0 1 

- . 300000000 3 788E+0 1 

. 207E-10 

9 

-.  2000000000057E+0 1 

-. 2000000000036E+0 1 

. 208E-10 

10 

-. 1 000000000004E+0 1 

-.9999999999878E+00 

. 160E-10 

1 1 

-.3882333782009E-I 1 

. 1472361028863E-10 

. 186E-10 

12 

. 9999999999943E+00 

. 1000000000003E+0 1 

.855E-1 I 

13 

. 2000000000046E+0 1 

. 2000000000053E+0 I 

. 726E- 1 1 

14 

. 3000000003798E+0 1 

. 3000000003806E+01 

. 800E-1 1 

15 

.4000000204992E+0 1 

. 400000020506 3E+01 

.710E-10 

16 

.50000081 5865 7E+01 

. 5000008 158664E+01 

.648E-1 1 

17 

.6000225680 1 67E+0 1 

.6000225680 1 87E+0 1 

.202E-10 

18 

. 7003952002645E+0 1 

.7003952002653E+01 

.821E-1 1 

19 

.8038941 1 19280E+01 

.80389411 19290E+0 1 

. 105E-10 

20 

. 92106 78647299E+01 

.92 106 7864 7308E+01 

. 875E- J 1 

21 

. 10746! 94 18286E+01 

,10746 19418288E+02  . 

. 143E-10 

Table  Vic:  m=300 

Eigenvalue  intervals  determined  by  EVSCAN  after  m Lanczos  iterations 
applied  to  W . 


I owe  i bounil  lal 


npporbonnd  (b) 


b-a 


1 


. iooioooooooookhu 

. 100 IOOOOOOOOOKH) 1 

. 782K- 1 1 

. IOI‘>OOOOOOOOOK.HU 

. 10  l‘»OOOOOOOOOK. HU 

. 924E- 1 3 

./OOOOOOOOOOOOKHlI 

. .’OOOOOOOOOOOOKHlI 

.48  IK- 12 

. 10  l‘)'l‘)‘)')‘)‘)')‘)(iK  f0  1 

. io  bi')')')>)‘)‘)«p»8i>o i 

. 1 J4K-  1 1 

Table  Vila: 

m“40  ops-l.  1*40*2 

see  see t i on  4 . 1 

lowotbounil  (a3 

uppe rbouinl  (b) 

b-rt 

. IOOIOOOOOOOOOKHU 

. I0(U000000000K>01 

. I42R-I3 

. 1 OO.’OOOOOOOOOK.  Ml  1 

. ioo/oooooooookhu 

./I  IK-14 

. ioo  looooooooo!-:  »o  i 

. 100  MlOOOOOOOOKHU 

. 42bK-  1 3 

. 1 004000000000K.H1 1 

. I004000000000KHU 

. 2 1 IK-  1 3 

. iooioooooooookhu 

. IOOIOOOOOOOOOKHU 

. 49 7K- 1 3 

. IOOpOOOOOOOOOKh) 1 

. IOObOOOOOOOOOK+OI 

. 426K-I  3 

. ioo/oooooooookhu 

. 100 /OOOOOOOOOK.HU 

. 81 3K- 1 3 

. IOOOOOOOOOOOOKhI 1 

. 1008000000000K.H1 1 

.426K.-1  3 

. IOO‘>OOOOOOOOOK  Hi  1 

. HHMOOOOOOOOOK.t-OI 

.71  IK- 14 

. ioiooooooooook.hu 

. lOIOOOOOOOOOOK.+OI 

. 426K- 1 3 

. 11'  1 1 000000000  IvfO  1 

.101  ioooooooook.hu 

. 284K- 1 3 

.10  1 .’000000000K  M)  1 

. 101  .’OOOOOOOOOKHU 

. 284K- 1 3 

. 101  1000000000K.  HU 

.101  ioooooooook.hu 

.476K-I3 

. i o 1 40000000001':  hi  t 

. I0I4000000000KHM 

. 7 1 3E- 1 .1 

. ioimkkiooooookhu 

. 10  1 SOOOOOOOOOK+O 1 

. 35.SK-1 3 

. IOIoOOOOOOOOOK.HU 

. UUbOOOOOOOOOKHU 

. 4')  7K-  1 3 

.101  /OOOOOOOOOKhU 

.101  /OOOOOOOOOK.HU 

. 2 1 3K- 1 3 

. 1018000000000KHM 

. 1018000000000K.hu 

. 3SSK- 1 3 

. IOI'>OOOOOOOOOKhU 

. 101 ‘1000000000K.hu 

. I42E-13 

. >0  l 

. .’OOOOOOOOOOOOK.hu 

. b‘)bK- 1 2 

. io.’0‘)')‘)>)8')>)«)‘),: hi  i 

. 1071  OOOOOOOOOKHU 

. 142K-13 

. 10.’ .’OOOOOOOOOK.HU 

. 102  .’OOOOOOl/OOKHU 

. 284 K- 1 3 

. io.’.’iM')‘)')>)‘)o>n>o  i 

. 10.’  1000000000K+0 1 

. 284K- 1 3 

. 10  .’4000000000K.HU 

. 1074000000000K.hu 

. 284K- 1 3 

. 10.’41>‘)9‘>‘>‘>'>‘)‘>KHU 

. tO/4')'l'i'»‘>‘i‘l‘)‘»KH>  1 

. 284K- 1 3 

. lO.’bOOOOOOOOOK.HU 

. W.’bOOOOOOOOOK+Ol 

. 284E-1 3 

. 10  .’/OOOOOOOOOK.HU 

. 10.’ /OOOOOOOOOK.HU 

. 42bK- 1 3 

. 10  .’8000000000K.hu 

. 3078000000000KH) 1 

0. 

. 10/8‘>‘)‘>‘>‘>‘)‘)‘>»)K  +0  1 

. 10.’8‘i‘i‘i')')‘)‘)')‘)KHU 

. I47K-I3 

. 10.’‘)‘>‘)‘)>)‘)ogi)«)KH)J 

. TO  10000000000KH)  1 

. 5b8K- 1 3 

. 10  10‘1‘1‘1‘l'WWK.HU 

. 10  30 ‘I 9 9 9 9 K.  HU 

. S68K-I3 

. 10  1 1 

. 10  n‘)‘)‘)‘»‘»)‘)‘)‘)K+01 

. 478F.-1 3 

. 10  1 f 0 1 

. to  l.,'l<)«)‘)>M>)0'JKH)  I 

. 284K-I3 

. 10  1 l>>*t'>‘>‘)‘W‘>‘>K.H)  t 

. 10  1 1909999999KH)  1 

.426K-13 

. 10  1.,‘)>)'»')')>)')>»)K  HU 

. 10T.',‘)‘l‘)‘>‘)‘i‘)‘)<)K.HU 

. I42E-I3 

. 10  1 

. 10  1S9‘19999999KH)  1 

.71  IK-13 

. 10  K>>»g*>‘»>>*)»)»)<iK.MU 

. 10  1<,990999999KHU 

.784K-I3 

. 10  1 /<)')')>)>|i)'l')')K  H)  1 

. 10  17990<I99999KH1 1 

. I47K.-I  3 

. 10  18‘)‘»‘»‘>‘»‘»‘>‘)9KHU 

. 10  lO'l'l'l'i'l'l'l'l'iK.fOl 

. 142E-1 3 

. 10  I'l'l ')*)'»'>')') iJOK.fl!  1 

. 1040000000000K.UU 

. I47K-I3 

Table  VI  lb: 

m“4  S eps  "4 , 7 *40 *2  ' ^ 

see  sect  ion  4.  1 
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U|'l>f  I tuMIUtl 


1 


i 


-.  ISO  1 lt>‘)7f /').,  U>0  I ISO  1 1 to)  7c, .’*>4  ll-X)  l 

I 7S7S')70  107  ISF.Ml  l -.1  7S/S')/0  107  I SIX'  I 

r.ililo  Villa:  , 7,0  l,anr<rt>*  .tto|>?*. 

IT  time  i o,|ui  i ,',t  tor  ISVl.AN:  IS.’  .tcviiJs 
I'l’-l  imo  i 1 fil  lor  K.VSllAN:  0 .*»  scMoiuln 


loWO  1 l>OIIII,l 

1 l'oitn.1 

. iso  it (.').•(. />>.,  it: mi  i 

-.11 

10  1 It. '!.’(• /')/,  IK.  M)  1 

. i /s  ’v»/o  ui;  i si: mi  i 

1 

1 S 7 S')  70  10.’  IMS. *0  1 

. 1 /-.  I')t>  It,.'  lOOfK.Ml  1 

1 

t')(,  u,.1  nun, Km'  » 

. 1 ’ IS‘>(.  ,.',it.  '.’M'.mI  I 

. 1 

' IS'lt, 77'>(>/7SK.M'  l 

. i / is  / s ; ts.'f.’K.Mi  i 

1 

• IS  /S/7  IS.'ti  .’K  M'  l 

. I / 17  , , ,7i>  ts/01  X)  1 

1 

' 1.’ is  70K  Ml  l 

. i / ii  i 07  so//  it- mi  i 

’ll.'  HI  .'Stl  7 7 IK  Ml  l 

. i //  is  it.:  Mso.,1  ui  i 

is  It,.'  /OS'l.’iK  Ml  1 

. 1 s ■ * ' /,,',  ■ ' Cl -* K Ml  1 

> ’ 7 7(>S  ’ 7t>/,t«K  M'  1 

. 1 '■(«')  7t>  ' t> S ’ '« 71'  Ml  1 

,(>')  (>  t>S  ‘ '< 7S.  Mi  1 

lalilr  VI  lllu  i-S|')*7  , t>0  Uim,.,o:i  Ht«'l>M. 

Ol’-I  into  r«'i|ui  t oil  tor  I SVl.AN:  .’.'ti  stu-ou,l:i 
01’— t imo  i i'ijiii  red  tor  KVSl’.ANt  0.*)  shrouds 


l OWi*  l I, omul 

i tnmrul 

. ISO  llt.'l  IK  Ml  1 

. ISO  llo'i.’.’sosOKM)  1 

. 1 IS 920  UI.'ISSMl  l 

. 1 7s  ",').’0  7'lS7,7K.Ml  l 

. 1 t‘lf>  l(,.’  lOOoKMl  1 

-.  1 /•',  i»c,  It,  7 lOOtiK  Ml  1 

. 1 / IS'lt,  .'SlMl  1 

-.17  IS'lt, 7 7'H, 7 7 SK. Ml  l 

. 1 / IS /S/7  IS.’A.’K.Ml  1 

-.17  IS7S77  IS7t,’7K. *0  1 

. 1 7 17,7, 7, 7f>  is  701  Ml  1 

1 7 !//,-,/, 7ti  IS70KM1  1 

-.17  117  10 7 SO  7 7 IK  M'  » 

.17  117  10  7 SO  7 7 IK.  Ml  1 

-.1/7  IS  It,  7 70S‘).’,KMl  1 

.1/7  is  It,  7 70 S') 7, K Ml  1 

. 1 7 701  It,  1 M 1 70SK  Ml  1 

. 1 7 70  1 It,  l')0*)/t,K.Ml  1 

. 1 7I-#  If,  1 /.#■>  IS.’KMl  1 

.17  1.,  If  1 ,S'»)  ISKM!  1 

.171  Ih-'iSOOSI  7 IK. Ml  1 

.171  It-’iM  7SOS7KM)  1 

. 1 S 7 7 '» S S -*  177  UK.MI  1 

. is  7 /'>S.,  ,t.S/,f,  IK.  Ml  1 

-.  IS  / 701  lt>0  7 1 77K  Ml  1 

.IS  / •OO'IS'l  M ASK  Ml  l 

. IS/S  70.,77'1S  17*  K M'  1 

. 1 s /',  '0  IS'l'i  77SK  Ml  l 

. IS/.,  1 lt,‘).,t,'l‘>'IKMl  l 

. 1 s/.,  1 It, '17, 000  IK  Ml  1 

. 1 S 7 7 7 7t«S  / 7h.,*>K  *0  1 

. 1 S 77  7 7t,S  7 7t,7,t,l  M'  l 

• . 1 st,')  7t,  7t,S  7 s7  7K  Ml  l 

. 1 St,')  7 1,  7t,S  7 S 7 /KM'  ' 

TaM,'  VI  lie: 

-SI  *>*10  ' # 60  l ,*nu* / o?4  Htepn 
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Gone  r a 1 inf o rma  1 1 on : 


- order  of  the  matrix  n-519. 

- the  matrix  was  available  on  a diskfile,  no  attempts  have  been 
made  to  optimize  the  matrix-vector  multiplications  required 
for  LSVLAN . 


- 60  - 


m 

1 

number  of 
eigenvalues 

eps 

992* 2"4 7 

P ” 

m 

number  of 
eigenvalues 

i 

eps 

992*  2~4  7 

30 

6 

6.7 

220 

56 

16. 

40 

8 

2.7 

230 

59 

16. 

50 

10 

2.9 

240 

62 

19. 

60 

14 

2.0 

250 

64 

13. 

70 

20 

2.8 

260 

65 

15. 

80 

23 

7.4 

270 

67 

17. 

90 

26 

6.2 

280 

68 

15. 

100 

28 

3.4 

290 

68 

18. 

1 10 

30 

6.4 

300 

73 

20. 

120 

32 

5.0 

310 

77 

28. 

130 

36 

11. 

320 

77 

11. 

140 

37 

8.4 

330 

77  ') 

116.  ') 

150 

38 

6.9 

340 

78 

19. 

160 

41 

9.4 

350 

83 

19. 

170 

46 

12.5 

360 

85 

10. 

180 

46 

7.6 

370 

87 

14. 

190 

49 

4.9 

380 

87 

20. 

200 

56 

16. 

390 

88 

17. 

210 

56 

13. 

400 

O' 

co 

9b.  ') 

Table  IX 

Results  of  EVSCAN  after  m I.anczog-itcrat  ions  applied  to  the  992-th 
order  matrix  described  rn  section  4.5 


’)  Results  after  second  scan  with  5*eps 


lowerbound 


upperbound 


-.23836 /2694904E+-00 

-. 238367269490 1E+00 

. 2383672694900F.+00 

. 2383672694904E+00 

Table  Xa: 

m- 1 5 

lowerbound 

upperbound 

-. 2383672694904K+00 

- . 2 383672694904E+00 

-.  14769007 I6874K+00 

-. 1 4769007 16874E+00 

-. I373857686379E+00 

-. I373857686379E+00 

. 1373857686379E+O0 

. 1373857686 379E+00 

. 14769007 168 74 E+00 

. 1476900716874E+00 

. 23836726949O3E+00 

. 2383672694903E+00 

Table  Xb: 

m=30 

lowerbound 

upperbound 

2383672694903K+00 

-. 2383672694903E+00 

-. 14769007 I6874E+00 

-. I4769007I6874E+00 

-. 1373857686379E+00 

-. 1373857686378E+00 

-. 10557 1057252 1E>00 

-. 10557 1057252 IE+00 

-. 101 7762399231E+00 

-. 101 776239923 1 E+00 

-.809883237301 2E-01 

-.80988323730 1 2E-0 1 

-.80148754 1 2532E-0 1 

-.80148754 1 2531E-0 1 

-.7608597598440E-01 

-. 7608597598440E-0 1 

. 7608597598434E-0 1 

. 7608597598435E-0 1 

.801487541 2526E-0 1 

.80148754 1 2526E-0 1 

. 8098832373009E-0 1 

.8098832373009E-01 

. 101 7762 399230E+00 

. I0I776239923OE+OO 

. 10557 105 72 520E+00 

. 10557 10572520E+00 

. 1373857686378E+00 

. I373857686378E+00 

. 14769007 1687 3E+00 

. 14769007 1 687  3E+00 

. 2383672694902E+00 

. 2383672694904E+00 

Tabic  Xc : m»60 

Eigenvalue  intervals  delivered  by  EVSCAN  after  m Lanczos- 
iterations  applied  to  the  product  matrix  described  in  4.6. 


Note:  The  values  in  the  table  above  should  be  multiplied  by 
i(-SQRT(-l))  so  that  they  represent  eigenvalues  of  A. 
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